Debugging: Changing ode equation during integration

Hello everyone,

I'm imprementing an adjustment to my chemical reactor design regarding physical limitation, due to the kinetics the value of one specie goes negative (comuptationally correct, but not physically) to account for the "no more fuel no more reaction" scenario i introduced a if:

ReactionRate(2*ovr) = ReactionK(2)*C(4)* 1e6 /CatalystDensity;  % kmol/kgcat/h
if y(ovr) < 0
    ReactionRate(2*ovr) = 0; % kmol/kgcat/h
end
dydt(i) = NetRateProduction(i)*(1-VoidFraction)*CatalystDensity*mw(i)/G;

in this way the dydt=0 and the value after reaching 0 (or a small negligible negative value) should stay constant, but after integration i still see a trend that goes way past 0

-0.0978571895684889

-0.0978554762967198

-0.0978537439506636

-0.0978519924638089

since the Ode for this specie depends only on the ReactionRate if dydt+0 it should stay constant, or is an error due to the fact that to evaluate if the rate should be 0 or not i use y? i tried with the last value of y using the persistent variable

persistent y_prev
.
if y_prev(ovr) < 0
    ReactionRate(2*ovr) = 0; % kmol/kgcat/h
end
%after the system of ODE
y_prev(ovr) = y(ovr)

how this happens? i tested with de-bugging to see if the Ode=0 if the condition on y is satisfied and eas zero, so a constant behaviour, but in the results the diminishing trend is still present.

Thank you so much for your help

[Torsten]

Try the "NonNegative" option of the ODE integrator instead of your construct from above.

And using more stringent tolerances ("RelTol" and "AbsTol") might help, too.

[Cris LaPierre]

It would be helpful if you could share a working example so we can reproduce the issue.