  The Liouville equation governing the evolution of the density matrix for an
atomic/molecular system is expressed in terms of a commutator between the
density matrix and the Hamiltonian, along with terms that account for decay and
redistribution. For finding solutions of this equation, it is convenient first
to reformulate the Liouville equation by defining a vector corresponding to the
elements of the density operator, and determining the corresponding
time-evolution matrix. For a system of N energy levels, the size of the
evolution matrix is N2xN2. When N is very large, evaluating the elements of
these matrices becomes very cumbersome. We describe a novel algorithm that can
produce the evolution matrix in an automated fashion for an arbitrary value of
N. As a non-trivial example, we apply this algorithm to a fifteen-level atomic
system used for producing optically controlled polarization rotation. We also
point out how such a code can be extended for use in an atomic system with
arbitrary number of energy levels.
