In short, I am interested in the simulation of large systems of nontrivial stochastic differential equations as they arise when applying phasespace methods to interacting quantum systems.
As an intro example: when one uses the socalled “positive P representation” on the twosite BoseHubbard model, this would result in eight nonlinear equations for the deterministic part, along with four independent Wiener processes that then enter the system as multiplicative noises, i.e noise_rate_prototype=zeros(8, 4)
.
Note: I am not asking how to solve this specific system.
Rather, I would like to have some input on the more general questions:

What algorithms are suggested for systems that have both, nonlinearity and multiplicative noise? (For the twosite BoseHubbard, I have tested
ImplicitEM()
, with satisfactory results in ‘unproblematic’ parameter regimes.) 
Numerically, how large can I make such systems and still get a reliable solution in finite time? For the combination of BoseHubbard and positive P, for instance, the degrees of freedom will have a factor of four w.r.t. the number of sites, and I will have half as many independent Wiener processes. Say I have 128 sites, that will make 512 degrees of freedom and 256 Wieners, combined in a coupled system of equations of the type described above. My naive intuition suggests that, among other things, one will have to create ever larger ensembles when the number of noises goes up. The nonlinearity will add further complications, obviously. Hence the issue of speed vs. accuracy seems to be of great concern here, too.
Thanks in advance for any help!