If I understand it correctly, I think the following is what you want: Solve the differential equation of `n` variables, with initial conditions defined using the previous differential equation solution of `n-1` variables and with an initial condition for the last variable (which randomly depends on the previous variable). The following code uses `iniCond` which generates initial conditions for a given `i`, `diffGen` generates equations with the initial conditions, and finally `solveDiffGen` solves the differential equations. iniCond[1] = {Subscript[x, 1][0] == 0.7}; diffGen[i_] := Join[Table[{Subscript[x, j]'[t] == Subscript[x, j][ t] (1 - Subscript[x, j][t] - nu (Sum[Subscript[x, k][t] Boole[k != j], {k, i}]))}, {j, i}], iniCond[i]] // Flatten; iniCond[i_] := Module[{conditions}, conditions = NDSolve[diffGen[i - 1], vars[[i - 1]], {t, 0, T}] /. (x_ -> g_) :> {(x /. {t -> 0}) -> (g /. {t -> T})} // Flatten; {conditions /. Rule -> Equal, {Subscript[x, i][0] == 0.01 Subscript[x, RandomInteger[{1, i - 1}]][t] /. t -> T}} // Flatten ]; solveDiffGen[i_] := NDSolve[diffGen[i], vars[[i]], {t, 0, T}];