  The density matrix is a widely used tool in quantum mechanics. In order to
determine its evolution with respect to time, the Liouville-von Neumann
equation must be solved. However, analytic solutions of this differential
equation exist only for simple cases. Additionally, if the equation is coupled
to Maxwell's equations to model light-matter interaction, the resulting
equation set -- the Maxwell-Bloch or Maxwell-Liouville-von Neumann (MLN)
equations -- becomes nonlinear. In these advanced cases, numerical methods are
required. Since the density matrix has certain mathematical properties, the
numerical methods applied should be designed to preserve those properties. We
establish the criterion that only methods that have a completely positive trace
preserving (CPTP) update map can be used in long-term simulations.
Subsequently, we assess the three most widely used methods -- the matrix
exponential (ME) method, the Runge-Kutta (RK) method, and the
predictor-corrector (PC) approach -- whether they provide this feature, and
demonstrate that only the update step of the matrix exponential method is a
CPTP map.
