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　non-associated　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.　©　2017,　The　Author(s).