Abstract: The loading/unloading judgment and analytical derivative operations have been the bottlenecks restricting numerical application of elastoplastic models. An unconstrained implicit stress update algorithm is proposed based on the hyper-dual step derivative approximation, which solves the above calculation difficulties. For the problem of loading/unloading judgment, in the new algorithm, the nonlinear stress integral equations with inequality constraints are transformed into an unconstrained minimization problem by using the smooth function to replace the Karush-Kuhn-Tucker conditions. Thus, there is no need for loading/unloading judgement during the calculation. To solve the problem of derivative evaluation, the algorithm uses the hyper-dual step derivative approximation instead of the analytical derivative to obtain the 1st derivative of the smooth function and the 1st and 2nd derivatives of the plastic potential function, which are used to construct iterative formulas for nonlinear calculation, ensuring the quadratic convergence speed of local stress update iterations and global equilibrium iterations. Numerical examples demonstrate that, compared with other numerical differentiation methods, the hyper-dual step derivative approximation is free from truncation errors and subtraction cancellation errors, and its computational results are almost equivalent to the analytical derivations. Finally, based on the proposed algorithm, a UMAT subroutine of smooth Mohr-Coulomb plasticity model is programmed. The effectiveness and convergence speed are verified through numerical analyses of three typical boundary value problems.