Simulating the nonlinear behavior of complex systems requires significant computational effort. Despite the rapid progress in computing technology, the demand is still strong for more efficient simulation methods in diverse structural dynamics fields such as nonlinear system identification and nonlinear system reliability. In addition to efficiency, algorithmic stability and accuracy must be addressed in the development of new simulation procedures. In this paper, a method to treat localized nonlinearities in a structure efficiently and accurately is proposed. The method is conditionally stable. The system equations are separated into a state-invariant linear part and a state-dependent nonlinear part that is considered to be external pseudo-forces that act on the linear system. The response of the system is obtained by fixed point iterations in which the pseudo-forces are updated until convergence. In addition to the one time-step-ahead prediction form, the novel idea of multiple time-step-ahead prediction is proposed. The effect of this approach is investigated and shown to increase algorithm efficiency and stability. To perform the numerical integration, time-stepping schemes like the exponential first-order hold method can be used to the advantage of efficiency and accuracy. To increase the accuracy and stability of the method, a novel second-order hold equivalent is derived and implemented. The efficiency, stability, and accuracy of the method are demonstrated in numerical examples. Finally, the method is applied to the earthquake-induced motion of a 20-story building with local nonlinearities. (C) 2016 Elsevier Ltd. All rights reserved.