# n-dimensional recursion

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.