Abstract: The equations of Magneto-hydrodynamics govern the flow of quasi-neutral electrically conducting fluids in the presence of induced and applied magnetic fields.  They are essentially the conservation equations from fluid mechanics strongly coupled with the Maxwell’s equations for the evolution of electric and magnetic fields.  The solution to these equations describe plasma phenomenon occurring in astrophysics, nuclear fusion, electric propulsion devices and in hypersonic flow control. The numerical development of a detailed Magneto-hydrodynamics (MHD) model will be presented. A matrix free implicit method is developed to solve the conservation equations within the framework of an unstructured grid finite volume formulation. The analytic form of the convective flux Jacobian is derived on a general unstructured mesh and used in a Lower-Upper Symmetric Gauss Seidel (LU-SGS) technique developed as part of the implicit scheme. A grid coloring technique is also developed to create data parallelism in the algorithm. The computational efficiency of the matrix free method is compared with two other methods: a global matrix solve technique that uses the GMRES (Generalized minimum residual) algorithm and an explicit method. The matrix-free method is observed to be overall computationally faster than the global matrix solve method and demonstrates excellent parallel scaling on multiple cores. The numerical model developed is subsequently used to study coaxial plasma thrusters and a comparison with experimental work is also presented.