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.