In this paper, we propose a new scheme for the numerical integration of the Landau–Lifschitz–Gilbert (LLG) equations in their full complexity, in particular including stray-field interactions. The scheme is consistent up to order 2 (in time), and unconditionally stable. It combines a linear inner iteration with a non-linear renormalization stage for which a rigorous proof of convergence of the numerical solution toward a weak solution is given, when both space and time stepsizes tend to \(0\). A numerical implementation of the scheme shows its performance on physically relevant test cases. We point out that to the knowledge of the authors this is the first finite element scheme for the LLG equations which enjoys such theoretical and practical properties.