Abstract-In this work, numerical analysis of nonlinear ferromagnetic problems is presented using the three-dimensional time-domain finite element method (TDFEM). Formulated with the secondorder nonlinear partial differential equation (PDE) combined with the inverse Jiles-Atherton (J-A) vector hysteresis model, the nonlinear problems are solved in the time domain with the NewtonRaphson method. To solve the ordinary differential equation (ODE) representing the magnetic hysteresis accurately and efficiently, several ODE solvers are specifically designed and investigated. To improve the computational efficiency of the Newton-Raphson method, the multi-dimensional secant methods, aka Broyden's methods, are incorporated in the nonlinear TDFEM solver. A nonuniform time-stepping scheme is also developed using the weighted residual approach to remove the requirement of a uniform time-step size during the simulation. The capability and the performance of the proposed methods are demonstrated by various numerical examples.