I am trying to set up a recursive function that generates `n`

number of differential equations for `Subscript[y, n][t]`

This function *almost* works.

`Table[{Subscript[y0, j] = 1}, {j, 50}];(*initial conditions for Subscript[y, n] assuming n<=50*) vars[n_] := {x, Table[Subscript[y, j], {j, n}]}; sol[0][T_, b_, d_, r_, n_] := sol[0][T, b, d, r, n] = Flatten[vars[n]] /. NDSolve[Flatten@{Join[{x'[t] == 10^7 r - d x[t] - b*x[t]*(Sum[Subscript[y, k][t], {k, n}])}, Table[Subscript[y, j]'[t] == -d Subscript[y, j][t] + b x[t]*Subscript[y, j][t], {j, n}], Flatten[Join[{x[0] == 0}, Flatten[Table[{Subscript[y, j][0] == Subscript[y0, j]}, {j, n}]]]]]}, Flatten[vars[n]], {t, 0, T}][[1]] sol[i][T_, b_, d_, r_, n_] := sol[i][T, b, d, r, n] = Flatten[vars[n]] /. NDSolve[Flatten@{Join[{x'[t] == 10^7 r - d x[t] - b*x[t]*(Sum[Subscript[y, k][t], {k, n}])}, Table[Subscript[y, j]'[t] == -d Subscript[y, j][t] + b x[t]*Subscript[y, j][t], {j, n}], Flatten[Join[{x[0] == 0},(*next bit seems to be the problem*) Flatten[Table[{Subscript[y, j][0] == sol[i - 1][T, b, d, r, n][[j + 1]][T]}, {j, n}]]]]]}, Flatten[vars[n]], {t, 0, T}][[1]] `

The initial condition `sol[0][T, b, d, r, n]`

works as expected and returns the interpolating functions:

`T = 4; b = 10^-7; d = 0.25; r = 0.2; n = 4; sol[0][T, b, d, r, n] `

And also returns the value at `t = T`

e.g. `sol[0][T, b, d, r, n][[2]][4]`

`sol[i][T, b, d, r, n] `

does not work for `i>0`

.

It seems that the problem is where the solution from the previous iteration `i-1`

is used to set the initial conditions for current iterate `i`

, as marked in the code.

I imagine this will be trivial to troubleshoot for someone on here. Any advice is much appreciated.