  We propose a novel algorithm for the temporal integration of the resistive
magnetohydrodynamics (MHD) equations. The approach is based on exponential
Rosenbrock schemes in combination with Leja interpolation. It naturally
preserves Gauss's law for magnetism and is unencumbered by the stability
constraints observed for explicit methods. Remarkable progress has been
achieved in designing exponential integrators and computing the required matrix
functions efficiently. However, employing them in MHD simulations of realistic
physical scenarios requires a matrix-free implementation. We show how an
efficient algorithm based on Leja interpolation that only uses the right-hand
side of the differential equation (i.e. matrix free), can be constructed. We
further demonstrate that it outperforms Krylov-based exponential integrators as
well as explicit and implicit methods using test models of magnetic
reconnection and the Kelvin--Helmholtz instability. Furthermore, an adaptive
step-size strategy that gives excellent and predictable performance,
particularly in the lenient- to intermediate-tolerance regime that is often of
importance in practical applications, is employed.
