In this paper the numerical approximation of stochastic differential equations satisfying a global monotonicity condition is studied. The strong rate of convergence with respect to the mean square norm is determined to be $$\frac{1}{2}$$ 1 2 for the two-step BDF-Maruyama scheme and for the backward Euler–Maruyama method. In particular, this is the first paper which proves a strong convergence rate for a multi-step method applied to equations with possibly superlinearly growing drift and diffusion coefficient functions. We also present numerical experiments for the $$\tfrac{3}{2}$$ 3 2 -volatility model from finance and a two dimensional problem related to Galerkin approximation of SPDE, which verify our results in practice and indicate that the BDF2-Maruyama method offers advantages over Euler-type methods if the stochastic differential equation is stiff or driven by a noise with small intensity.