You can generate correlated gaussian variables
by computing the Cholesky decomposition of the variance matrix
and multiplying independent gaussian variables with it.
There was a similar question on Quant Stackexchange
a few days ago.
To simulate the solution of a stochastic differential equation (SDE),
you can just discretize it, i.e., consider it as a difference equation (this is called the "Euler scheme") and use the correlated
brownian motion you want.
In some cases, this is too imprecise,
but it is possible to derive
a second-order approximation of the SDE
from Ito's formula (this is called the "Milstein scheme";
it is an analogue of the Runge-Kutta method).
If you are using R, this is implemented (for instance)
in the sde
package
and described in the accompanying book.