A robust and concise implicit stress integration algorithm of elastoplastic models is presented. It does not require the loading/unloading estimation and analytical derivation operation for the stress update. First, the elastoplastic stress update problem is recast into an unconstrained minimization problem by utilizing the smooth function to bypass the loading/unloading estimation. Then, the object problem is solved by the line search method instead of the Newton method for better convergence. The consistent tangent operator is evaluated by the complex step derivative approximation without the subtraction cancellation error, which provides the quadratic convergence rate of global iteration. The rationality of the numerical consistent tangent operator is validated by the one obtained by the analytical derivation. A recently developed non‐orthogonal elastoplastic (NEP) clay model is implemented using the new algorithm. The algorithm is confirmed through comparing the numerical solution and the analytical one for a cavity expansion problem. The algorithm performance is assessed based on a series of geotechnical boundary value problems. It is found that the new algorithm is more robust than the one employed by ABAQUS. The source code of the model implementation can be downloaded from https://github.com/zhouxin615.