We show how to simulate a model of many molecules with both strong coupling to many vibrational modes and collective coupling to single-photon modes. This is done by combining the process tensor matrix product operator method with a mean-field approximation that reduces the dimensionality of the problem. We analyze the steady state of the model under incoherent pumping to determine the dependence of the polariton lasing threshold on cavity detuning, light-matter coupling strength, and environmental temperature. In addition, by measuring the correlation twice, we investigate the second order variation with respect to the mean field and calculate the photoluminescence spectrum. Our method allows us to simulate many-body systems with strong coupling to multiple environments and extract both static and dynamic properties.
