We propose a new robust, accurate, and fast numerical method for solving the Landau-Lifshitz equation which describes the relaxation process of the magnetization distribution in ferromagnetic material. The proposed numerical method is second-order accurate in both space and time. The approach uses the nonlinear multigrid method for handling the nonlinearities at each time step. We perform numerical experiments to show the efficiency and accuracy of the new algorithm on two- and three-dimensional space. The numerical results show excellent agreements with exact analytical solutions, the second-order accuracy in both space and time, and the energy conservation or dissipation property. (C) 2014 Elsevier B.V. All rights reserved.