This paper describes a robust and efficient elasto-plastic boundary element method (BEM) intended for complex problems in geomechanics. Some novel techniques are employed to evaluate the singular integral equations accurately and efficiently. In particular, the strong singularities arising in the domain integrals of the stress integral equations are removed by transforming them from domain integrals to (cell) boundary integrals. Adaptive integration schemes are used throughout to optimize the integration process. The system equations are condensed by recasting them in terms of the plastic multiplier, and a novel Newton–Raphson algorithm delivers greatly improved solution times. Some examples demonstrate the scope and power of the method.