Hi everyone,

I have a recurring problem in my ODE-based models. In most case the outputs of the ODE system have to be strictly positive, but because the precision of the solver is not perfect, outputs that should be 0 are sometimes estimated as small negative values. This leads to problems when I use these outputs in further computations, e.g. ratios may become negative or greater than 1, sqrt transformations don’t work, etc.

Until now I used a silly way around, checking for negative values and replacing them by a small positive value:

```
y = integrate_ode_rk45(ode, init, t0, theta, xrs, xis, 1.0E-7, 1.0E-7, 1.0E3);
for(i in 1:S) for(j in 1:C) y[i,j] = y[i,j]<0 ? 1e-10 : y[i,j];
```

I found that forcing y to be positive leads to errors and interrupts sampling.

Is there a more principled way to deal with this issue?