We developed a new theoretical approach to simulate polariton linear and non-linear spectra, which fully account for the non-Markovian exciton-phonon couplings and the effect of cavity loss. This approach is based on combining the partial linearized density matrix (PLDM) approach, which is a non-adiabatic path-integral formalism, and a new Lindblad method that exactly maps the Lindblad dynamics with a stochastic evolution of the mapping variables. We refer to this approach as the L -PLDM method. Using this approach, we simulated the linear and third-order response functions, and upon Fourier transform, obtain the corresponding optical spectra.