  We describe a generalized algorithm for evaluating the steady-state solution
of the density matrix equation of motion, for the pump-probe scheme, when two
fields oscillating at different frequencies couple the same set of atomic
transitions involving an arbitrary number of energy levels, to an arbitrary
order of the harmonics of the pump-probe frequency difference. We developed a
numerical approach and a symbolic approach for this algorithm. We have verified
that both approaches yield the same result for all cases studied, but require
different computation time. The results are further validated by comparing them
with the analytical solution of a two-level system to first order. We have also
used both models to produce results up to the third order in the pump-probe
frequency difference, for two-, three- and four-level systems. In addition, we
have used this model to determine accurately, for the first time, the gain
profile for a self-pumped Raman laser, for a system involving 16 Zeeman
sublevels in the D1 manifold of 87Rb atoms. We have also used this model to
determine the behavior of a single-pumped superluminal laser. In many
situations involving the applications of multiple laser fields to atoms with
many energy levels, one often makes the approximation that each field couples
only one transition, because of the difficulty encountered in accounting for
the effect of another field coupling the same transition but with a large
detuning. The use of the algorithm presented here would eliminate the need for
making such approximations, thus improving the accuracy of numerical
calculations for such schemes.
