This paper develops a steady-state simulation algorithm for magnetic devices using the shooting-Newton method. We show how any diagonally-implicit Runge-Kutta method may be used for numerical integration. Calculation of the Jacobian of the nonlinear state transition function is avoided by employing GMRES to calculate the correction to the solution. An interior permanent magnet synchronous machine is simulated to examine the effects of higher order numerical integration methods on the accuracy of the electromagnetic torque, stator winding losses, and rotor permanent magnet losses. For the problem under consideration, a great decrease in simulation time is achieved for a given accuracy.