  We present a multi-dimensional numerical code to solve isothermal
magnetohydrodynamic (IMHD) equations for use in modeling astrophysical flows.
First, we have built a one-dimensional code which is based on an explicit
finite-difference method on an Eulerian grid, called the total variation
diminishing (TVD) scheme. Recipes for building the one-dimensional IMHD code,
including the normalized right and left eigenvectors of the IMHD Jacobian
matrix, are presented. Then, we have extended the one-dimensional code to a
multi-dimensional IMHD code through a Strang-type dimensional splitting. In the
multi-dimensional code, an explicit cleaning step has been included to
eliminate non-zero $\nabla\cdot B$ at every time step. To estimate the
proformance of the code, one- and two-dimensional IMHD shock tube tests, and
the decay test of a two-dimensional Alfv\'{e}n wave have been done. As an
example of astrophysical applications, we have simulated the nonlinear
evolution of the two-dimensional Parker instability under a uniform gravity.
