# How to replace the expression in nest?

I need to calculate the following function with replacement in recurrence

f[x0_, y0_]:= (s1 + s2)/t1 /. {s1 ->
NDSolveValue[{x''[t] + x[t] == 0., y''[t] + x[t]^2 y[t] == 0.,
x[0.] == x0, x'[0.] == 0., y[0.] == y0, y'[0.] == 0.},
y, {t, 0, t1}][t1], s2 -> NDSolveValue[{x''[t] + x[t] == 0, y''[t] + x[t]^2 y[t] == 0.,
x[0.] == x0, x'[0.] == 0., y[0.] == y0, y'[0.] == 1.},
y, {t, 0, t1}][t1]} /. {t1 -> Take[Reap[
NDSolve[{x''[t] + x[t] == 0., x[0.] == x0, x'[0.] == 0.,
WhenEvent[x'[t] > 0., {Sow[t], "StopIntegration"}]},
x, {t, 0., 100.},
MaxStepSize -> 0.001]], {2, -1}][[1]][[1]][[1]]}


But I meets an error says "NDSolveValue: Endpoint t1 in {t, 0., t1} is not a real number".

The key expression is

NDSolveValue[{x''[t] + x[t] == 0, y''[t] + x[t]^2 y[t] == 0.,
x[0.] == 1., x'[0.] == 0., y[0.] == 1., y'[0.] == 0.},
y, {t, 0, Re[t1]}][Re[t1]] /. {t1 -> Take[Reap[
NDSolve[{x''[t] + x[t] == 0, x[0.] == 1., x'[0.] == 0.,
WhenEvent[x'[t] > 0., {Sow[t], "StopIntegration"}]},
x, {t, 0., 100.},
MaxStepSize -> 0.001]], {2, -1}][[1]][[1]][[1]]}


which comes the error. How can I solve this problem?

• It is unclear to me what you are trying to achieve. The innermost equation has an analytical solution ($x(t)=\cos{t}$), for instance. Commented Dec 11, 2020 at 17:02
• You are asking Mathematica to numerically solve for s1 with a symbolic endpoint t1 and then afterwards replacing t1 with a calculated value. To avoid the warning, compute t1 first. Commented Dec 11, 2020 at 19:18
• @SimonWoods, because I construct a function with this expression, so I need it compute in one time. Commented Dec 12, 2020 at 0:18
• @MarcoB, I construct a function which is similar to the expression. The simple example here is just to show you the error I meet, not the exactly expression I need. Commented Dec 12, 2020 at 0:27
• When I execute your code, I get no error. Please include the code that causes the error. Commented Dec 12, 2020 at 1:08

Make the t1 replacement in the s1,s2 substitution, grouping with parentheses:

ff[x0_, y0_] := s1 + s2 /. (
{s1 :>
NDSolveValue[{x''[t] + x[t] == 0., y''[t] + x[t]^2 y[t] == 0.,
x[0.] == x0, x'[0.] == 0., y[0.] == y0, y'[0.] == 0.},
y, {t, 0, t1}][t1],
s2 :> NDSolveValue[{x''[t] + x[t] == 0,
y''[t] + x[t]^2 y[t] == 0., x[0.] == x0, x'[0.] == 0.,
y[0.] == y0, y'[0.] == 1.}, y, {t, 0, t1}][t1]} /. {t1 ->
NDSolveValue[{x''[t] + x[t] == 0., x[0.] == x0, x'[0.] == 0.,
WhenEvent[x'[t] > 0., {"StopIntegration"}]},
Indexed[x["Domain"], {1, -1}], {t, 0., 100.},
MaxStepSize -> 0.001]}
);

ff[-1, -3]

(*  -2.54137  *)


Alternative:

ff[x0_, y0_] :=
Block[{t1 =
NDSolveValue[{x''[t] + x[t] == 0., x[0.] == x0, x'[0.] == 0.,
WhenEvent[x'[t] > 0., {"StopIntegration"}]},
Indexed[x["Domain"], {1, -1}], {t, 0., 100.},
MaxStepSize -> 0.001]},
s1 + s2 /.
{s1 :>
NDSolveValue[{x''[t] + x[t] == 0., y''[t] + x[t]^2 y[t] == 0.,
x[0.] == x0, x'[0.] == 0., y[0.] == y0, y'[0.] == 0.},
y, {t, 0, t1}][t1],
s2 :> NDSolveValue[{x''[t] + x[t] == 0,
y''[t] + x[t]^2 y[t] == 0., x[0.] == x0, x'[0.] == 0.,
y[0.] == y0, y'[0.] == 1.}, y, {t, 0, t1}][t1]}
];

• Thank you so much! Commented Dec 12, 2020 at 1:30
• What's the function of the extra "Indexed[x["Domain"], {1, -1}]" in your expression? Commented Dec 12, 2020 at 1:35
• And it doesn't calculate correctly when I change the definition of the function into ff[x0_, y0_] := (s1 + s2)/t1 Commented Dec 12, 2020 at 1:41
• @JieJiang Use Block instead of replacement: f[x0_, y0_] := Block[{t1 = NDSolveValue[...]}, ...]. You can look up Indexed. It's a delayed version of Part. Or maybe you mean x["Domain"], which gives the interval of the domain of an InterpolatingFunction. Commented Dec 12, 2020 at 1:57
• See mathematica.stackexchange.com/questions/28337/… for more about InterpolatingFunction, which is what NDSolve computes. Commented Dec 12, 2020 at 1:58