  We describe a numerical scheme for studying time-dependent, multifluid,
magnetohydrodynamic shock waves in weakly ionized interstellar clouds and
cores. Shocks are modeled as propagating perpendicular to the magnetic field
and consist of a neutral molecular fluid plus a fluid of ions and electrons.
The scheme is based on operator splitting, wherein time integration of the
governing equations is split into separate parts. In one part independent
homogeneous Riemann problems for the two fluids are solved using Godunov's
method. In the other equations containing the source terms for transfer of
mass, momentum, and energy between the fluids are integrated using standard
numerical techniques. We show that, for the frequent case where the thermal
pressures of the ions and electrons are << magnetic pressure, the Riemann
problems for the neutral and ion-electron fluids have a similar mathematical
structure which facilitates numerical coding. Implementation of the scheme is
discussed and several benchmark tests confirming its accuracy are presented,
including (i) MHD wave packets ranging over orders of magnitude in length and
time scales; (ii) early evolution of mulitfluid shocks caused by two colliding
clouds; and (iii) a multifluid shock with mass transfer between the fluids by
cosmic-ray ionization and ion-electron recombination, demonstrating the effect
of ion mass loading on magnetic precursors of MHD shocks. An exact solution to
a MHD Riemann problem forming the basis for an approximate numerical solver
used in the homogeneous part of our scheme is presented, along with derivations
of the analytic benchmark solutions and tests showing the convergence of the
numerical algorithm.
