Methods for modeling large driven-dissipative quantum systems are becoming increasingly urgent due to recent experimental progress in a number of photonic platforms. We demonstrate the positive-P method to be ideal for this purpose across a wide range of parameters, focusing on the archetypal driven-dissipative Bose-Hubbard model. Notably, these parameters include intermediate regimes where interactions and dissipation are comparable, and especially cases with low occupations for which common semiclassical approximations can break down. The presence of dissipation can alleviate instabilities in the method that are known to occur for closed systems, allowing the simulation of dynamics up to and including the steady state. Throughout the parameter space of the model, we determine the magnitude of dissipation that is sufficient to make the method useful and stable, finding its region of applicability to be complementary to that of the truncated Wigner method. We then demonstrate its use in a number of examples with nontrivial quantum correlations, including a demonstration of solving the urgent open problem of large and highly nonuniform systems with tens of thousands of sites
