Nested NDSolve with WhenEvent: setting up new equations and discarding old ones

Question

I am checking a complex situation where a function x1[t1] only exit in a finite time span, and after this a new function x2[t2] will be set up according to the final value of x1, and then repeat the same procedure to build x3[t3], ..., xn[tn].

A most simple 3-nest toy code can be rawly formatted like this

d0 = 1;
{c1, c2, c3} = {1, 2, 3};
end=20;

res1 = NDSolve[{
   x1'[t1] == c1, x1[0] == 0,
   WhenEvent[x1[t1] == d0,
    res2 = NDSolve[{
       x2'[t2] == c2, x2[t1] == x1[t1],
       WhenEvent[x2[t2] == 2*d0,
        res3 = NDSolve[{
           x3'[t3] == c3, x3[t2] == x2[t2]
           }, x3, {t3, t2, end}, MaxSteps -> 1000000];
        "StopIntegration"
        ]}, x2, {t2, t1, end}, MaxSteps -> 1000000];
    "StopIntegration"
    ]}, x1, {t1, 0, end}, MaxSteps -> 1000000];

We can now check res1, res2 and res3, and they work. But I hope that this code can be more concise than its current state, like how we discard the previous functions and code the remaining inner-nest part, especially for multiple-nest cases, ie. a 100-nest case. I think the most proper way would be to go meta-programming together with recursion... am I right?

---

Update

Based on Albert Retey's answer, we can furthe treat x[t] as a vector function, like this

cVals = {{1, 1/2}, {1/3, 1/4}, {1/5, 1/10}};
tvals = {5, 10, 20};
tstart = 0;
xsol = Quiet[
   NDSolveValue[{
     x'[t] == cVals[[n[t]]], x[tstart] == {0, 0}, n[tstart] == 1,
     WhenEvent[t == tvals[[n[t]]],
      If[TrueQ[n[t] < Length[cVals]], n[t] -> n[t] + 1, "RemoveEvent"]]
     }, x, {t, tstart, 20},
    DiscreteVariables -> {Element[n, Integers]}, MaxSteps -> 1000000
    ], Part::pspec];

Plot[xsol[t], {t, ##}] &[Sequence @@ xsol["Domain"][[1]]]

Where we used cVals = {{1, 1/2}, {1/3, 1/4}, {1/5, 1/10}}.

I think the code can be more general if we can use something like cVals = {1, {1/2, 1/3}, {1/4, 1/5, 1/6}} so that the x[t] can have variable dimensions at different time scopes. I wonder if this idea is feasible?

---

Update 2

Here is a first try, which shows the idea, but the code does not work yet...

cVals = {1, {1/2, 1/3}};
tvals = {5, 10};
tstart = 0;
xsol = Quiet[
   NDSolveValue[{
     x'[t] == cVals[[n[t]]], x[tstart] == 0, n[tstart] == 1,
     WhenEvent[t > tvals[[n[t]]],
      If[TrueQ[n[t] < Length[cVals]],
       xv = x[t];
       n[t] -> n[t] + 1;
       x[t] -> {xv, xv},
       "RemoveEvent"]]
     }, x, {t, tstart, 20},
    DiscreteVariables -> {Element[n, Integers]}, MaxSteps -> 1000000
    ], Part::pspec];

Plot[xsol[t], {t, ##}] &[Sequence @@ xsol["Domain"][[1]]]

Answer 1

As mentioned in my comment, it is possible for the example you have shown to achieve the same thing without a recursion and without even restarting NDSolve. The trick is to introduce a discrete variable which is changed at each event. Here is something that I think does the same thing as your code:

cVals = {1, 2, 3};
tstart = 0;
xsol = Quiet[NDSolveValue[{
     x'[t] == cVals[[n[t]]], x[tstart] == 0, n[tstart] == 1,
     WhenEvent[x[t] == n[t]*1,
      If[TrueQ[n[t] < Length[cVals]], n[t] -> n[t] + 1, "RemoveEvent"]
      ]
     },
    x, {t, tstart, 20},
    DiscreteVariables -> {Element[n, Integers]},
    MaxSteps -> 1000000
    ], Part::pspec];

Plot[xsol[t], {t, (xsol@"Domain")[[1, 1]], (xsol@"Domain")[[1, 2]]}]

the code should work for lists of arbitrary length for cVals. The Quiet for Part::pspec is necessary because the equations will be evaluated for symbolic values of n[t]. It could be avoided by defining e.g. cc[n_?NumericQ]:=cVals[[n]] and use that in the first equation.

For other cases you might not be able to reformulate the switch in such a way, but then the first thing I'd try would be a loop, here is code which does that:

c = {1, 2, 3}; d0 = 1; tstart = 0; tend = 20;

solution[n_, nmax_] := Module[{t},
   x[n] = NDSolveValue[{
      x[n]'[t] == c[[n]], x[n][tstop[n - 1]] == x[n - 1][tstop[n - 1]],
      If[n == nmax,
       Unevaluated[Sequence[]],
       WhenEvent[x[n][t] == n*d0, "StopIntegration"]
       ]
      },
     x[n], {t, tstop[n - 1], tend}, MaxSteps -> 1000000
     ];
   tstop[n] = (x[n]@"Domain")[[1, 2]];
   ];

tstop[0] = tstart; x[0][tstart] = 0;

Do[solution[n, Length[c]], {n, 1, Length[c]}]

This checks that the solutions are really calculated as desired:

DownValues[x] // TableForm

Compared to writing the same thing by nesting this won't run into any iteration limit problems. Anyway, for completeness here is code which does what I think is closest to what was asked for originially, a recursive function:

c = {1, 2, 3};
d0 = 1;
x[0][0] = 0;
tend = 20;

solution[n_, nmax_, tprev_] := Module[{t},
   res[n] = NDSolve[{
       x[n]'[t] == c[[n]], x[n][tprev] == x[n - 1][tprev],
       If[n == nmax,
        Unevaluated[Sequence[]],
        WhenEvent[x[n][t] == n*d0, solution[n + 1, nmax, t]; 
         "StopIntegration"]
        ]
       },
      x[n], {t, tprev, tend}, MaxSteps -> 1000000
      ];
   ];

solution[1, Length[c], 0];

DownValues[res]