A stress update algorithm for elastoplastic models with nonconvex yield surfaces is presented. An explicit algorithm based on the Runge-Kutta embedded method of second-order accuracy is developed. The crossing of the yield surface is properly taken in to account by means of a robust intersection-finding algorithm. This algorithm is based on a simple multiple-root-finding technique, which requires the yield function, evaluated along a given secant stress path, to be continuously differentiable to the second order. The accuracy of the intersection-finding algorithm is illustrated through examples using simple yield functions similar to the ones adopted in models for unsaturated soils and modified Cam clay model yield surface with nonconvex Argyris et al. failure criterion. Isoerror surfaces and finite element simulations are used to investigate the accuracy and efficiency of the stress update algorithm. It is observed that although the algorithm for nonconvex surfaces is slower than its equivalent for convex surfaces, the accuracy can be controlled locally by means of specified tolerances. Copyright (c) 2008 John Wiley & Sons, Ltd.