# Undetected WhenEvent with NDSolve

I was considering the following piece of code

u[y_, z_] = 36*y*(1 - y)*z*(1 - z);
sys = {x'[t] == u[y[t], z[t]], y'[t] == 0, z'[t] == -2.2*^-3,
x[0] == 0, y[0] == y0, z[0] == z0,
WhenEvent[x[t] > 500, "StopIntegration"]};
sol = Table[
NDSolve[sys, {x, y, z}, {t, 0, 500}], {y0, 0.1, 0.9, 0.1}, {z0,
0.1, 0.9, 0.1}];
ParametricPlot3D[{x[t], y[t], z[t]} /. sol, {t, 0, 500},
PlotRange -> {{0, 1000}, {0, 1}, {0, 1}}, ImageSize -> 500,
BoxRatios -> {5, 1, 1}]


However, it does not seem that WhenEvent captures correctly the event x>500, indeed Mathematica does not stop integration. I've also tried with "==" but I have the same problem.

What is wrong?

Thanks!

• I think it's actually working. You're just plotting the solutions for t from 0 to 500 even for the solutions that stop at some earlier time. Try looking at just the solution for {y0 -> 0.5 , z0 -> 0.9} and do x["Domain"]/.sol to get the domain of the function. Also, you can replace "StopIntegration" with Print[{y0, z0, t}]; "StopIntegration" to make Mathematica spit out the times at which it stops and for which parameters. Oct 22, 2015 at 16:02
• The core issue here is Plot (et al) suppress the warnings about extrapolation of the InterploatingFunction. Is there a way to disable that (ie. show the warnings )? Oct 22, 2015 at 16:56
• @March: Yes, I've realized just now that it is a problem of extrapolation outside the integration interval. However, Print[{y0,z0,t}] does not print anything, furthermore x["Domain"]/.sol gives {0,500} and it seems to me that the integration does not stop at all. SORRY! I was unaware that when you put in more events these must be listed! Actually it is working, thanks! Oct 22, 2015 at 18:38
• Any clue how to Reap/Sow in WhenEvent? Oct 22, 2015 at 18:45

u[y_, z_] := 36*y*(1 - y)*z*(1 - z);
sys = {
x'[t] == u[y[t], z[t]], y'[t] == 0,
z'[t] == -2.2*^-3,
x[0] == 0, y[0] == y0, z[0] == z0,
WhenEvent[x[t] == 500, "StopIntegration"]};

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];
sol = Table[NDSolveValue[sys, {x, y, z}, {t, 0, 500}], {y0, 0.1, 0.9, 0.2},
{z0, 0.1, 0.9, 0.2}];
doms = InterpolatingFunctionDomain[#[[1]]][[1, 2]] & /@  Flatten[sol, 1];
tt = Transpose[{Flatten[sol, 1], doms}];
plots = ParametricPlot3D[Through[#[[1]][t]], {t, 0, #[[2]]},
BoxRatios -> {5, 1, 1}] & /@ tt;
prs = PlotRange /. AbsoluteOptions[#, PlotRange] & /@ plots;
prs1 = {Min[#[[1]]], Max[#[[2]]]} & /@ Transpose /@ Transpose@prs;
Show[plots, PlotRange -> prs1]


edit

You may plot the ending time wrt y0 and z0:

coords = Table[{y0, z0}, {y0, 0.1, 0.9, 0.1}, {z0, 0.1, 0.9, 0.1}];
prs = (Last@First@(PlotRange /. AbsoluteOptions[#, PlotRange])) & /@  plots;
ListPlot3D[Join @@@ Transpose@{Flatten[coords, 1], List /@ prs}]


• Interesting solution though complicate. I have to study in detail, thanks. Oct 22, 2015 at 18:35