systems of ODE with positive solution constraints

Hello everybody!

I'm here to ask help to the MHF community regarding a problem about systems of ODE I'm trying to solve manually (C++). As it appears to be nightmarish, I'm hoping there exists some standard solution that I'm unable to find (either C++ or Matlab libraries, or in the worst case - but still possibly better than mine - some standard algorithm I can try to implement).

I'm dealing with a chemical problem, in which many species are present at some concentration (variables x1, x2, ...) and some reactions are ongoing (differential equations in the form x1'=f1(x1,x2,...), x2'=f2(x1,x2,...), ...). I want to be able to sketch the kinetic evolution of the system in the time frame [0,t0]. The problem I'm facing is that *the concentrations cannot go negative*. Hence when the Eulero method calculates a decrease in the variable, say, x1 which is greater than the value of x1 itself, I need x1 to just go to zero instead. Also, if you have at the same time x2'=x1 and x3'=2x1 in two separate differential equations, and x1 is almost zero, x1 needs to go to zero while at the same time increasing x3 twice as much as x2. This would be all easy if it were not for the fact that if you consider having x3'=x1 and x3'=x1+x2, and x1 is going to disappear in the next step, then x1 must split equally in the two pathways (becoming x3 alone, and becoming x3 together with x2) ONLY if there is enough x2; if there isn't, then it will no longer be split equally, it will get consumed by x2 as much as possible, hence becoming x3 on two pathways (x1->x3 and x1+x2->x3), but not any longer in a 1:1 ratio. (Wondering)

In other words, I have a system of ODEs and the solutions should not become negative. How can I implement a numerical solution that retains the different ratios of production of different species, with the solution not becoming negative, AND considering the possibilities that some variables may be limiting and may therefore be changing the ratios at which changes of variables will be operated in every time step?

I'm sorry for explaining things a bit between chemistry and mathematics and I hope that, if you have ideas to share, you will ask me the extra details that are needed, since I will explain them, hopefully more precisely than I already did! (Happy)

Thanks in advance!

Elio