A finite-difference time-domain formalism for simulating coherent linear pulse propagation is presented that incorporates a medium response described by any two-time energy gap correlation function. Two algorithms, for real and complex correlation functions, are developed to evaluate the electric polarization through explicit treatment of the density matrix for a two-level system. The coherence relaxation terms in the resulting finite-differenced Maxwell-Liouville equations depend on integrals over the energy gap fluctuation correlation function. The algorithms are used to simulate ultrashort mid-infrared pulse propagation through optically dense samples of HDO in liquid D2O as a demonstration of their performance and flexibility. These algorithms represent a first step toward the goal of incorporating complicated material responses into the full-field simulation of nonlinear pulse propagation and nonlinear optical spectroscopy. (C) 2003 American Institute of Physics.