We present an accurate, efficient, and flexible method for propagating spatially distributed density matrices in anharmonic potentials interacting with solvent and strong fields. The method is based on the Nakajima-Zwanzig projection operator formalism with a correlated reference state of the bath that takes memory effects and initial/final correlations to second order in the system-bath interaction into account. A key feature of the method proposed is a special parametrization of the bath spectral density leading to a set of coupled equations for primary and N auxiliary density matrices. These coupled master equations can be solved numerically by representing the density operator in eigenrepresentation or on a coordinate space grid, using the Fourier method to calculate the action of the kinetic and potential energy operators, and a combination of split operator and Cayley implicit method to compute the time evolution. The key advantages of the method are: (1) The system potential may consist of any number of electronic states, either bound or dissociative. (2) The cost for arbitrarily long solvent memories is equal to only N + 1 times that of propagating a Markovian density matrix. (3) The method can treat explicitly time-dependent system Hamiltonians nonperturbatively, making the method applicable to strong field spectroscopy, photodissociation, and coherent control in a solvent surrounding. (4) The method is not restricted to special forms of system-bath interactions. Choosing as an illustrative example the asymmetric two-level system, we compare our numerical results with full path-integral results and we show the importance of initial correlations and the effects of strong fields onto the relaxation. Contrary to a Markovian theory, our method incorporates memory effects, correlations in the initial and final state, and effects of strong fields onto the relaxation; and is yet much more effective than path integral calculations. It is thus well-suited to study chemical systems interacting with femtosecond short laser pulses, where the conditions for a Markovian theory are often violated.