  We present a new numerical model for solving the Chew-Goldberger-Low system
of equations describing a bi-Maxwellian plasma in a magnetic field.
Heliospheric and geospace environments are often observed to be in an
anisotropic state with distinctly different parallel and perpendicular pressure
components. The CGL system represents the simplest leading order correction to
the common isotropic MHD model that still allows to incorporate the latter's
most desirable features. However, the CGL system presents several numerical
challenges: the system is not in conservation form, the source terms are stiff,
and unlike MHD it is prone to a loss of hyperbolicity if the parallel and
perpendicular pressures become too different. The usual cure is to bring the
parallel and perpendicular pressures closer to one another; but that has
usually been done in an ad hoc manner. We present a physics-informed method of
pressure relaxation based on the idea of pitch-angle scattering that keeps the
numerical system hyperbolic and naturally leads to zero anisotropy in the limit
of very large plasma beta. Numerical codes based on the CGL equations can,
therefore, be made to function robustly for any magnetic field strength,
including the limit where the magnetic field approaches zero. The capabilities
of our new algorithm are demonstrated using several stringent test problems
that provide a comparison of the CGL equations in the weakly and strongly
collisional limits. This includes a test problem that mimics interaction of a
shock with a magnetospheric environment in 2D.
