|Abstract The hyperbolic function proposed by Abbo- Sloan was employed not only to approach the Mohr- Coulomb criterion but also to express the plastic potential function. A better approximation to the Mohr-Coulomb yield and potential surfaces was achieved by increasing the transition angle and proven to be highly efficient in numerical convergence. When a Gaussian integral point goes into plastic state, two cases on yield stress adjustments were introduced. They may avoid solving the second derivative of the plastic potential function and the inverse matrix compared with the existing subroutine. Based on the above approaches, a fully implicit backward Euler integral regression algorithm was adopted. The two- and three-dimensional user subroutines which can consider the associated or non-associated flow rule were developed on the platform of the finite element program—ABAQUS. To verify the reliability of these two subroutines, firstly, the numerical simulations of the indoor conventional triaxial compression and uniaxial tensile tests were performed, and their results were compared with those of the embedded Mohr-Coulomb model and the analytical approach. Then the main influential factors including the associated or nonassociated flow rule, the judgment criteria of slope failure, and the tensile strength of soil were analyzed, and the application of the two-dimensional subroutine in the stability analysis of a typical soil slope was discussed in detail through comparisons with the embedded model and the limit analysis method, which shows that this subroutine is more applicable and reliable than the latter two.