I've recently got a question from a Master student about a numerical simulation/integration of the SDE of the following shape $$ \mathrm dX_i(t) = \left(\sum_{j=1}^M \nu_{ji} a_jX_j(t)\right)\mathrm dt + \sum_{j=1}^M\left(\nu_{ji}\sqrt{a_j X_j(t)}\right)\mathrm dW_j(t) $$ where $\nu_{ji}$ are entries of the integer-valued matrix (stoichiometry matrix), and $a_j$ are positive reals (rates of chemical reactions). The student tried the Euler-Maruyama scheme $$ X_i(t+\Delta t) = X_i(t) + \left(\sum_{j=1}^M \nu_{ji} a_jX_j(t)\right)\Delta t + \sum_{j=1}^M\left(\nu_{ji}\sqrt{a_j X_j(t)\Delta t}\right)\xi_j $$ where $\xi_j$ are iid standard normal random variables. However, for $\Delta t$ small enough the scheme does not converge and blows up drastically. This is of course expected: for small $\Delta t$ the term in front of $\xi_j$ significantly dominates the deterministic drift term.
We also looked into higher-order schemes for the integration, such as Runge-Kutta method, but they all blow up due to similar reasons. I would be happy to hear of any hint how to approach this problem.