A numerical method for solving the Landau–Lifshitz–Gilbert equation for an ensemble of superparamagnetic nanoparticles within the mean-field approximation is presented. The classical fourth-order Runge–Kutta method is employed for the time integration of the equation. The model simulates an ensemble of uniaxial nanoparticles subjected to a constant external magnetic field. It is shown that the proposed approach accurately reproduces the magnetization dynamics: the components perpendicular to the field decay, while the longitudinal component relaxes toward a steady-state value. The results are qualitatively consistent with previously published data obtained using the Vinamax simulation software.