Implicit calculation of elastoplastic models involves the construction of the Jacobian matrix of nonlinear stress integral equations in local stress update iterations and the consistent tangent operator in global equilibrium iterations, both of which require a substantial number of 1st and 2nd derivatives. For simple constitutive models, the derivative terms can be obtained by the manual differentiation method which is accurate and efficient. However, the manual differentiation method can quickly become difficult, time-consuming, and error-prone as the constitutive models increase in complexity. The numerical differentiation methods are increasingly favored by researchers since such methods are easy to implement and not sensitive to the complexity of the functions. This chapter 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 plasticity 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.