A two-dimensional axisymmetric computational algorithm is developed to simulate the plasma flow field in a self-field MPD thruster, in order to determine the flow behavior and the electromagnetic characteristics distribution. The convective flux vector is computed by using Roe's scheme in combination with Powell's eigensystem technique, and a new modified MUSCL technique called OMUSCL2 is employed to obtain the stable high-accuracy solution. Madrane-Tadmor entropy correction is used to prevent unrealistic expansion shocks near the electrodes tips. To accurately capture the physics of plasma in the system, different physical-chemical sub-models including multi-level non-equilibrium ionization model, generalized Ohm's law for partially ionized plasma, micro-instabilities effects, two-temperature model, and a real equation of state are considered. Numerical results of plasma flow simulation in a cylindrical lab-scale thruster, with mass flow rate of 6 g/s and total discharge current of 8 kA, are presented and comparison with experimental data shows good agreement between the predicted and measured contours of enclosed current and electric potential. The estimated thrust is 16.34 N which exhibits less than 5% difference compared with measured value. Furthermore, this simulation properly predicts the experimentally observed argon jet structure. (C) 2014 IAA. Published by Elsevier Ltd. All rights reserved.