We show how to simulate a model of many molecules with both strong coupling to many vibrational modes and collective coupling to a single photon mode. We do this by combining process tensor matrix product operator methods with a mean-field approximation which reduces the dimension 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. Moreover, by measuring two-time correlations, we study quadratic fluctuations about the mean field to calculate the photoluminescence spectrum. Our method enables one to simulate many-body systems with strong coupling to multiple environments, and to extract both static and dynamical properties.