A geomechanics program for wellbore stability analysis has been developed consisting of two modules: an analytical-based solution and a numerical-based solution. In the first part, input data are imported, including petrophysical well logs, pressure data, formation well tops, and a well path. Lithology intervals are set with proper prediction equations to calculate rock mechanical properties based on laboratory tests. In-situ stress and pore pressure are determined using different methods, including the poroelastic plane strain model and stress polygon. From the theory of plane strain, new equations are solved to determine horizontal tectonic strains (ϵ h, ϵ H) from drilling events such as total mud loss and breakout during drilling. Safe mud weight bounds are calculated through depth and in different azimuths and inclinations applying the Mohr-Coulomb and the Mogi-Coulomb failure criteria. The latter underestimated the minimum mud weight to prevent wellbore breakout. The transversely vertical isotropy of shale formation is programmed with multiple stress transformations via the weak-plane method. In the second module, a 3D model around the wellbore is discretized with hexahedral eight-point elements and programmed using the finite-element (FE) method. Rock mechanical property and displacement boundary conditions are applied to solve FE equations. Stress from the numerical model matched to the Kirsch model and results show that maximum stress concentration around the wellbore corresponds to the wellbore breakout, which has analytically been established. A new well plan across the 3D model was examined to obtain the safe mud weight bounds and results were in agreement with the analytical calculations.