The positive-P method for bosons has been shown to allow for scalable exact simulations of many-body open quantum systems, provided dissipation is strong enough to stabilise stochastic trajectories. Designing similar methods for spin systems, relevant to quantum technologies, has been an ongoing challenge. Recently, we developed a novel Positive-P-based method for systems with
spin-1/2 degrees of freedom by mapping them onto sets of three stochastic complex variables. This approach combines seamlessly with the usual bosonic positive-P method and so is applicable to systems of coupled spins and bosons. Unlike previous attempts at positive-P representations based on SU(2) coherent states for spins, our formalism restores the important advantage of positive-P simulations in that cases that are stable converge to the exact quantum mechanical solution. In this talk, I will briefly introduce the established formalism for bosonic positive P, before explaining our new approach that allows us to treat spins in a related way. I will present results and analysis of the regimes of stability of the method for a driven dissipative Jaynes-Cummings model, and show its scalability to an example of an extended driven dissipative Jaynes-Cummings-Hubbard lattice.