Solvers for Differential Equations are not converging (Julia)

Question

Problem description

I am trying to get numerical solutions for the stochastic schrödinger equation. A simplified version of the equation that I am solving for the ostensible state vector x:

dx/dt = -i*A*alpha(t)*x

where A is an arbitrary matrix and alpha(t) is a complex Gaussian noise processes with zero mean and given covariance and arbitrary pseudocovariance (given the total kernel is positive semidefinite). For example the covariance of alpha(t) looks like:

E[alpha(t) alpha(t')] = e^{-i\omega(t-t')}

Methodology

To construct the noise process I generated the whole (timesteps x timesteps)-dim covariance matrix and did an eigenvector eigenvalue decomposition to get the matrix B so that B*B^{\dagger} = e^{-i\omega t}. Then I constructed the process alpha(t) = B * X with X a Gaussian variable with covariance as the identity matrix and mean 0. Then alpha(t)should have the desired properties. (Is there a better method to generate the desired colored noise? I also tried with alpha(t) = beta_1(t) + i*beta_2(t) and beta real valued processes, and the covariance E[(beta_1(t) + i*beta_2(t))(beta_1(t) - i*beta_2(t))] = cos(\omega(t-t')) - i*sin(\omega(t-t')), with the same problems that I will explain now.)

The physical property, the density operator, can be computed as average over the outer product of different trajectories of the ostensible state vector x. The trace of the density operator should stay constant at 1 (with some fluctuations due to the noise). I used different solvers for (R)ODE but for all of them the trace is diverging, for some of them very rapidly to ~10^20, for some to ~100, but for none of them it is staying constant as it should. I used 100,000 trajectories which should be enough.

For the (R)ODE I first generated a trajectory for the noise, then passed it to the correct shape and then passed it to the RODESolver:

my_noise = NoiseGrid(t, noise_trajectory)

# Solve the RODE
prob_RODE = RODEProblem(f, u0, tspan, p, noise = my_noise)
sol_RODE = solve(prob_RODE, RandomEM(), dt=dt)

For ODE I also computed the trajectory first but then took alpha(t) as time-varying parameter.

noise_trajectory = fun_noise_trajectory()  # computes a trajectory for the noise
    
function f(du, u, p, t)
    parameter_index = Int(floor(t / dt)) + 1  # Adjust index calculation
    my_noise = noise_trajectory[parameter_index]  # Get parameter value at the index
    du .= - 1im * A * u .* my_noise
end;

p = (noise_trajectory, )

# Solve the ODE
prob_ODE = ODEProblem(f, u0, tspan, p)
sol_ODE = solve(prob_ODE, RK4(), dt=dt)

Questions

I think the problem comes from the different solvers, since the solutions differ a lot for them. RODE works better than ODE when comparing with the analytical solutions.

What does RODE do differently than ODE, if one also passes a already computed trajectory to the solver?

What is the RandomEM() solver doing numerically (it is the only implemented solver for RODE?

Which method should I use to get an accurate result?