I have a system of ODEs for functions $f_i(t)$. At each time $t$, $f_i(t)$ is the position of particle $i$. The functions $f_i$ have a monotonicity property: at all times $0 < f_1(t) < f_2(t) < \dots < f_n(t)$. This property is preserved by the true dynamics, since there is a divergent repulsive term in the ODEs when any two $f_i$ get close together. But they do, in fact, get close together. In fact they get close enough together that crossings can easily happen in a numerical simulation (and if they do then the rest of the run is ruined).
I have been simulating my system with ode15s in MATLAB (all the other solvers fail quite spectacularly). My method has been to use an events function to search for collisions, and if I find a collision, then I just let the simulation fail. (This is reasonably convenient to program, because I need the events function anyway for an unrelated reason.) I am not sure how efficient this approach is. Would it be more efficient to construct ODEs for $g_i=f_i-f_{i-1}$ (with the convention $f_0=0$) and then utilize the NonNegative option, then just reconstruct $f_i$? Or would that end up being implemented in basically the same way internally? (I don't really know much about how odezero works.)