SUMMARYFor numerical simulations to be e ective tools in plasma propulsion research, a high-order accurate solver that captures MHD shocks monotonically and works reliably for strong magnetic ÿelds is needed. For this purpose, a characteristics-based scheme for the MHD equations, with ux limiters to improve spatial accuracy, has been developed. In this method, the symmetric form of the MHD equations, accounting for waves propagating in all directions, are solved. The required eigensystem of axisymmetric MHD equations, with appropriate normalization, is presented. This scheme was validated with unsteady (Riemann problem) and force-free equilibrium (Taylor state) test cases, as well as with measured current density patterns in a magnetoplasmadynamic thruster.