Numerically solving SDE in Python

Looking for reference on doing numerical solving of a system of SDE's in python. I am trying to model a stochastic SIR model:

$$\begin{cases} dS = -\beta SIdt - \sigma SIdW \\ dI = (\beta SI - \gamma I)dt + \sigma SIdW \\ dR = \gamma Idt \end{cases}$$

I have been using a python package called SDEint, however, I am running into errors of getting negative sample paths which should not occur.

I have posted my code here: Stochastic SIR using SDEint python package

I have thought of trying to build my own solver but I am not terrible proficient in python. Any references to examples or literature that will help me in constructing my own solver is appreciated.

Additionally, perhaps there is a way to parametrize the system to avoid getting negative values in the sample paths. Such a parametrization is mentioned here How is the simulation done of the black scholes model?

You can really ignore $R$ in this setting; if it is desired it can be reconstructed as $R=S_0+I_0+R_0-S-I$.

With that in mind, you can look at $f(S,I)=\ln(S)$ and $g(S,I)=\ln(I)$. You can use Ito's formula on each of these. For $f$ (assuming I didn't make any mistakes) you have

$$df=\left ( S^{-1} (-\beta S I) - \frac{1}{2} S^{-2} \sigma^2 S^2 I^2 \right ) dt - S^{-1} (\sigma S I) dW \\ =\left ( -\beta I - \frac{1}{2} \sigma^2 I^2 \right ) dt - \sigma I dW$$

and for $g$ you have

$$dg=\left ( I^{-1} (\beta S I) - I^{-1} (\gamma I) - \frac{1}{2} I^{-2} \sigma^2 S^2 I^2 \right ) dt + I^{-1} (\sigma S I) dW \\ =\left ( \beta S - \gamma - \frac{1}{2} \sigma^2 S^2\right ) dt + \sigma S dW$$

and now you can replace $S$ with $e^f$ and $I$ with $e^g$.

It is important here that $S$ and $I$ and thus $f$ and $g$ are driven by the same Brownian motion, so that the noise term conserves the total population.

With this formulation there is no more barrier to be concerned about crossing. The reconstructed solutions for $S$ and $I$ will not cross $0$ by definition since they will be written as $e^f$ and $e^g$ respectively.