I came across a system of differential equations in the form: $\newcommand{\D}[1]{\frac{\mathrm{d}#1}{\mathrm{d}x}}$ \begin{align} f_1(x,y,z)\D{y}+f_2(x,y,z)\D{z}&=f_3(x,y,z),\\ f_4(x,y,z)\D{y}+f_5(x,y,z)\D{z}&=f_6(x,y,z). \end{align}
I numerically solved the system as follows. Let $\Phi$ be given by $$ \Phi(x,y,z)=\begin{bmatrix}y'\\z'\end{bmatrix}= \begin{bmatrix} f_1(x,y,z)&f_2(x,y,z)\\ f_4(x,y,z)&f_5(x,y,z) \end{bmatrix}^{-1}\begin{bmatrix} f_3(x,y,z)\\ f_6(x,y,z) \end{bmatrix}, $$
then Euler's method takes the form $$ \begin{bmatrix} y\\ z \end{bmatrix}_{(i+1)}= \begin{bmatrix} y\\ z \end{bmatrix}_{(i)} +h~\Phi(x,y,z). $$
I tested this approach (using RK4 instead of Euler) for simple systems (for which I can obtain analytical solutions) and it seems to work.
Other than the matrix of $\{f_1,f_2,f_4,f_5\}$ being singular, are there any other conditions for which this method fails?