# Solving ODEs produces error NDSolve::nrnum1

I'm a new member but I have been reading the great answers here for sometime now. I am trying to solve a system of ODEs for various values of a parameter epsilon, using a slider. I am using WhenEvent to start one of the variables, Z, only once another variable, T, has reached a certain value.

It seems to work OK except that I get the following error:

NDSolve::nrnum1: "The function value -2000. + 1.50707*10^-9 i is not a real number when the arguments are {357.662, 6.12468*10^-7 + 1.50707*10^-9 i, -6.53218*10^-11 - 1.50707*10^-9 i, 8213.09 + 0. i}"

I can't figure out why I am getting this error. Any suggestions would be greatly appreciated.

My code is:

Clear[Y, T, Z, Gamma, Mu, k, r, Epsilon, a, m, Rho];

eqnY[Y_, Z_] := γ (1 - Y/k) Y - μ Y - r ε Sqrt[Y] - a Y Z;
eqnT[Y_, Z_, T_] := r ε Sqrt[Y] - a T Z;
eqnZ[Z_] := ρ Z - m Z;

system = {
Y'[t] == eqnY[Y[t], Z[t]],
T'[t] == eqnT[Y[t], T[t], Z[t]],
Z'[t] == eqnZ[Z[t]],
Y[0] == 1, T[0] == 0, Z[0] == 0
};
control = {WhenEvent[T[t] == 2000, Z[t] -> 1]};

params =
{k -> 8000, γ -> 68, μ -> 0.2, a -> 0.01, m -> 0.01, r -> 0.05, ρ -> 0.1, ε -> 0.001};

Manipulate[
Block[{tmax = 500, sol},
sol =
NDSolve[{system, control} /. {ε -> εin} /. params,
{Y, T, Z}, {t, 0, tmax},
MaxSteps -> 5^6];
Plot[Evaluate[{Y[t], T[t], Z[t]} /. sol], {t, 0, tmax},
PlotRange -> {{0, tmax}, {0, 10000}},
Frame -> {True, True, False, True},
Axes -> False,
FrameLabel -> {"Time", "Densities"},
PlotLegends -> {"Y", "T", "Z"},
ImageSize -> {500}]],
{{εin, 0.01, "ε"}, 0, 5, Appearance -> "Labeled"},
ContinuousAction -> False]

-
you're getting the error because, as it plainly says, there is a small imaginary part. You could try Chop –  acl Mar 21 at 13:39