  We present a method to perform a variational analysis of the quantum master
equation for driven-disspative bosonic fields with arbitrary large occupation
numbers. Our approach combines the P representation of the density matrix and
the variational principle for open quantum system. We benchmark the method by
comparing it to wave-function Monte-Carlo simulations and the solution of the
Maxwell-Bloch equation for the Jaynes-Cummings model. Furthermore, we study a
model describing Rydberg polaritons in a cavity field and introduce an
additional set of variational paramaters to describe correlations between
different modes.
