In quantum optics, it is often necessary to simulate the equations of motion of a system coupled to a reservoir. Using a Schrö dinger picture approach, this can be done either by integrating the master equation for the density matrix[1] or by