Accurately　evaluating　derivatives　poses　a　key　challenge　when　numerically　implementing　complex　constitutive　models.　This　work　presents　an　implicit　stress　update　algorithm　that　utilizes　the　hyper　dual　step　derivative　approximation　to　address　derivative　evaluations　in　elastoplastic　problems.　Initially,　the　performance　of　various　numerical　differentiation　methods　is　discussed　and　compared　by　examining　their　numerical　errors　in　the　representative　example.　Subsequently,　the　hyper　dual　step　derivative　approximation,　without　truncation　and　subtractive　cancellation　errors,　is　employed　to　compute　the　Jacobian　matrix　and　consistent　tangent　operator,　ensuring　quadratic　convergence　in　both　local　and　global　computations.　The　size　of　the　Newton　search　step　is　optimized　by　the　line　search　technique,　thereby　enhancing　the　convergence　in　solving　nonlinear　stress　integral　equations.　Finally,　the　proposed　stress　update　algorithm　is　used　to　implement　the　non-associated　Mohr–Coulomb　plastic　model　in　the　ABAQUS　software　using　the　UMAT　subroutine.　The　stress　update　algorithm＇s　performance　and　its　practical　application　in　geotechnical　engineering　problems　are　demonstrated　using　five　boundary　value　problems.　©　2023　Elsevier　B.V.