  Generalized master equations valid for the third order response of an
optically driven multi-level electronic system are derived within Zwanzig
projection formalism. Each of three time intervals of the response function is
found to require specific master equation and projection operator. Exact
cumulant response functions for the harmonic profiles of the potential energy
surfaces leading to Gaussian spectral diffusion are reproduced. The proposed
method accounts for the nonequilibrium state of bath at the border of
intervals, and can be used to improve calculations of ultrafast non-linear
spectra of energy transferring systems
