In this study, a two-dimensional, axisymmetric, computational Algorithm has been developed to simulate the plasma flowfield in a MPD thruster for the purpose of determining the flow behavior and electromagnetic characteristics distribution. The solution employs Roe’s flux vector difference method in combination with Powell’s characteristics-splitting scheme. To ensure the stable high-accuracy solution, new modification of MUSCL technique so called OMUSCL2 method is used. According to being supersonic strong gasdynamic expansion near the electrodes tip, HHT entropy correction is employed. Further improvements to the physical model, such as the inclusion of relevant classical transport properties, a real equation of state, multi-level equilibrium ionization models, anomalous transport, and multi-temperature effects, that are essential for the realistic simulation MPD flows, are implemented. Numerical results of a lab-scale thruster are presented, whereby comparison with experimental data shows good agreement between the predicted and measured enclosed current and electric potential.