I would like to simulate the following system of interacting OU processes on $[0,T]$:
$$dX_t^1=(X_t^2-X_t^1)\,dt+\sigma_1 \,dW_t^1,\quad X_0^1=x_1$$
$$dX_t^2=(X_t^1-X_t^2)\,dt+\sigma_2 \,dW_t^2,\quad X_0^2=x_2$$
where $W^1$ and $W^2$ are two independent Wiener processes. I am aware of the fact there is a closed-form formula for the classical one-dimensional OU process, but in the case of two interacting processes as above I don't believe this is possible. As a result, my first idea was to naively apply the Euler-Maruyama scheme, as follows:
$$X_{t_{k+1}}^1=X_{t_k}^1+(X_{t_k}^2-X_{t_k}^1)h+\sigma_1(W_{t_{k+1}}^1-W_{t_k}^1)$$
$$X_{t_{k+1}}^2=X_{t_k}^2+(X_{t_k}^1-X_{t_k}^2)h+\sigma_2(W_{t_{k+1}}^2-W_{t_k}^2)$$
with $X_{t_0}^1=x_1$, $X_{t_0}^2=x_2$, $t_k=\frac{kT}{N}$, N being the number of subdivisions, and $h=\frac{T}{N}$.
However, I am not sure whether the Euler-Maruyama scheme can be applied to this interacting case, and whether or not the scheme would effectively converge in this specific context. Any ideas or references to literature would be greatly appreciated, thanks.
You can simply try the transformation $$A_t = \frac{X_t^1 + X_t^2}{\sqrt{2}}, \quad B_t = \frac{X_t^1 - X_t^2}{\sqrt{2}}$$ The equations will get much simpler.