The MHD equations have served plasma physicists, space physicists and astrophysicists very well for the past several decades. Even so, one comes to the realization that the MHD equations are only an approximation that applies to systems where the collision rate between charged particles is fast enough to rapidly restore isotropy in the pressure tensor. For many systems involving collisionless plasmas, this assumption of the isotropy of the pressure tensor is invalid. In such situations, the Chew Goldberger and Low [1] equations take over. They are based on a double adiabatic approximation and allow the pressure parallel to the field lines to differ from the pressure perpendicular to the field lines. While the CGL equations are well-formed, they have resisted numerical solution for a very long time. The reason is that if these equations are taken by themselves, they could potentially give rise to an unbounded increase in pressure anisotropy in certain circumstances. This unbounded increase in anisotropy can eventually result in a loss of hyperbolicity and a consequent code crash. Real world systems, like the solar wind, or tokamaks, show that there are limits to the pressure anisotropy. Therefore, to make sense of the CGL equations, and their physical import, one has to understand the plasma physics effects, such as the mirror instability and firehose instability, which restore pressure isotropy at a microphysical level. Folding in these physical effects also results in numerical challenges that have to be overcome. In this talk we show how the CGL equations can be made numerically tractable by a deft combination of physical reasoning and advanced numerical methods. Preliminary applications will also be presented.