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 matrix where about three quarters of the elements vanish, 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 the presence of noise and inertia, and demonstrate an excellent agreement between simulations and theory for the phase space distribution.