System of SDE with complex noise

Let’s say I want to solve numerically the following stochastic partial differential equation:
\begin{equation} du(x,y,t)=\nabla^{2}u(x,y,t) dt+k\cdot dW^{r}(x,y,t)+i k\cdot dW^{c}(x,y,t)) \end{equation},
where u is a complex function in 2d, i is the imaginary unit, and k is a constant.
As one can see, we have a complex noise. In particular, the Winer process has real part W^{r} and imaginary part W^{c}, which are independent random variables.
Now, in order to solve the equation, I discretize in space, such that I obtain the following system of equations:


I didn’t explicitely write the drift part (the laplacian) to make it shorter, and also because I want to concentrate on the stochastic part, as I think it is the cause of my problems.
So let me rewrite just the stochastic part of the system of equations:

du[i,j]=k\cdot dWr[i,j]+i k \cdot dWc[i,j]
(the “i” multiplying k outside the squared brackets is the imaginary unit and not an index!).
Now, reading the page of DifferentialEquations.jl ()DifferentialEquations.jl documentation, I understood that julia handles the stochastic part of an SDE in the following form:
\begin{equation} du_{stochastic}=G*dW \end{equation}
where, if u is a vector, then G is a matrix such that the matrix-vector product G*dW gives the stochastic change of the solution du. (So I guess the solver multiplies G by an array of independent noises dW)

Now, in my case u is a matrix, so in order to put myself in the same scenario described in the documentation, I reshaped u from a matrix (N \times M) into a column vector (N M \times 1) ( usually I have N=M=300).
Having done that, the G matrix in my case will be N M \times 2N M, such that, if I multiply it by a vector of independent noises dW of dimension 2N M \times 1, I get a vector of dimension N M \times 1 corresponding to du. In particular, by looking at my system of equations, my G matrix will have the following form:
\begin{equation} g=\begin{bmatrix} k & i \cdot k & 0 & 0 & 0 & 0 & \dots \\ 0 & 0 & k & i \cdot k& 0 & 0 & \dots \\ 0 & 0 & 0 & 0 & k & i \cdot k & \dots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \\ \end{bmatrix} \end{equation}.

Therefore, I implemented the noise prototype as follows:

function create_matrix_G_prototype(imax,jmax,p)
           for i in 1:imax*jmax 
    return  sparse(g_matrix_prototype)


In this way, I can implement the G function to pass to the solver as an in-place function:

function G!(du, u, p, t)

Then I define the problem and solve it with:


The issue here is that, first of all, the code seems to be terribly slow. Moreover, I get a warning telling me that the dt chosen by the solver is smaller than the dt_{min}:

┌ Warning: dt(8.881784197001252e-16) <= dtmin(8.881784197001252e-16) at t=0.37906220534331775, and step error estimate = 1.1208004104491986. Aborting. There is either an error in your model specification or the true solution is unstable.

This is surprising to me since I didn’t expect the equation to be stiff (but at this point I’m starting to guess it is). Has anyone any idea of what is going wrong? what shall I do about the warning? maybe choose a specific solver for stiff equations?
Also, is there a better approach that I could take for this specific case in order to increase the performance and speed of the code? Is it a good idea to separate real and imaginary parts instead of dealing with complex values?