This paper studies the use of time stepping methods to integrate the equations which model the magnetodynamic fields in ferromagnetic laminations. A new high-order time stepping method is compared to well-known time stepping methods (Euler, Crank Nicolson, ...) for a diffusion problem of electromagnetic fields in magnetic laminations with hysteresis. For validation purposes, firstly a theoretical problem with analytical solution is solved. Afterwards, more realistic problems are studied. As the presented method is A-stable, accurate numerical loops can be computed, even in the case of hysteresis effects.