How to replace the expression in nest?

Question

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?

Answer 1

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]}
   ];