Radiative cooling plays a crucial role in the dynamics of many astrophysical flows, and is particularly important in the dense shocked gas within Herbig-Haro (HH) objects and stellar jets. Simulating cooling processes accurately is necessary to compare numerical simulations with existing and planned observations of HH objects, such as those from the Hubble Space Telescope and the James Webb Space Telescope. In this paper, we discuss a new, non-equilibrium cooling scheme we have implemented into the three-dimensional magnetohydrodynamic (MHD) code ASTROBEAR. The new cooling function includes ionization, recombination, and excitation of all the important atomic species that cool below 10 000 K. We tested the routine by comparing its predictions with those from the well-tested one-dimensional Cox–Raymond shock code (Raymond 1979). The results show that ASTROBEAR accurately tracks the ionization fraction, temperature, and other MHD variables for all low-velocity (≲90 km s−1) magnetized radiative shock waves. The new routine allows us to predict synthetic emission maps in all the bright forbidden and permitted lines observed in stellar jets, including H α, [N II], [O I], and [S II]. We present an example as to how these synthetic maps facilitate a direct comparison with narrowband images of HH objects.