A new approach is proposed for nonlinear boundary element methods in computing internal stresses accurately using a complex-variable formulation. In this approach, the internal stresses are obtained from the numerical derivatives of the displacement integral equations that involve only weakly singular integrals. The collocation points in the displacement integral equations are dened as complex variables whose imaginary part is a small step size for numerical derivatives. Unlike the finite difference method whose solution accuracy is step-size dependent, the complex-variable technique can provide ‘‘numerically-exact’’ derivatives of complicated functions, which is step-size independent in the small asymptotic limit. Mean while, it also circumvents the tedious analytical differentiation in the process. Consequently, the evaluation of the nonlinear stress increment only deals with kernels no more singular than that of the displacement increment. In addition, this technique can yield more accurate stresses for nodes that are near the boundary. Three examples are presented to demonstrate the robustness of this method.