  We present a novel method to compute the phase space distribution in the
nonequilibrium stationary state of a wide class of mean-field systems involving
rotators subject to quenched disordered external drive and dissipation. The
method involves a series expansion of the stationary distribution in inverse of
the damping coefficient; the expansion coefficients satisfy recursion relations
whose solution requires computing a sparse matrix, making numerical evaluation
simple and efficient. We illustrate our method for the paradigmatic Kuramoto
model of spontaneous collective synchronization and for its two mode
generalization, in presence of noise and inertia, and demonstrate an excellent
agreement between simulations and theory for the phase space distribution.
