A general numerical method is presented for the Liouville-von Neumann integro-differential equation of motion of a reduced density matrix p, for molecular systems which arise when delayed (non-Markovian) dissipative dynamics are considered. Our method is a fourth-order extended Runge-Kutta integration scheme, which can be generalized for use with matrices. The method is applied to a spin-boson model. A comparison is made with results in the limits of instantaneous dissipation and Markovian dissipation. (c) 2005 Elsevier B.V. All rights reserved.