  We present a novel method to solve the Maxwell-Liouville-von Neumann (MLN)
equations in an accurate and efficient way without invoking the rotating wave
approximation (RWA). The method is a combination of two established concepts,
namely the operator splitting method as well as the adjoint representation of
the Lie algebra SU(N) (or pseudospin representation). The former concept
ensures the accuracy of the approach, but is computationally expensive. The
latter concept provides an efficient representation of the problem and two
optimization possibilities. We have implemented and verified both optimization
approaches and demonstrate that substantial speedup can be achieved.
