How to perform integration processes NDSOlve and show list of random variables used in this process?

Question

I have an ODE system which solves 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 one of the previous variables r[n_]). In addition, the system has a other random variable (\[Xi][j_, k_] different of r[n_]). The code of the system is shown below:

nmax = 5;(*max number of equations*)
T = 70;(*Max Time*)
\[Xi][j_, k_] := RandomReal[];(*random variable*)
r[n_] := RandomInteger[{1, 
   n}];(*set up unknown vars and differential equations,for n \
species*)
vars := Table[Subscript[x, j], {j, n}];
eqns := Table[
   Subscript[x, j]'[t] == 
    Subscript[x, j][
      t] (1 - Subscript[x, j][
        t] - (Sum[\[Xi][j, k] Subscript[x, k][t] Boole[k != j], {k, 
          n}])), {j, n}];
ics = {Subscript[x, 1][0] == 0.7};(*initial ICs*)
(*main loop*)
Do[(*solve for n species*)
 sol[n] = NDSolve[{eqns, ics}, vars, {t, 0, T}][[1]];
 (*set up ICs for n=n+1 species*)
 ics = Join[
   Table[Subscript[x, j][0] == 
     Evaluate[Subscript[x, j][T] /. sol[n]], {j, 
     n}], {Subscript[x, n + 1][0] == 
     Evaluate[0.01 Subscript[x, r[n]][T] /. sol[n]]}];, {n, nmax}]

I try to do two things:

i) Do 1000 different process as described above;

I thought about adding a Do as:

Do[Do[(*solve for n species*)
 sol[n] = NDSolve[{eqns, ics}, vars, {t, 0, T}][[1]];
 (*set up ICs for n=n+1 species*)
 ics = Join[
   Table[Subscript[x, j][0] == 
     Evaluate[Subscript[x, j][T] /. sol[n]], {j, 
     n}], {Subscript[x, n + 1][0] == 
     Evaluate[0.01 Subscript[x, r[n]][T] /. sol[n]]}];, {n, nmax}],{pp,1000}]

But I was not successful.

ii) Obtain a list with all the values of the variable \[Xi][j, k] drawn in the 1000 operations.

Can someone help me?

Answer 1

1. Your attempt in (i) fails because you forgot to initialize ics for every round of repeation.
2. Generating the list of ξ[j, k] before the summation will make the collection of it easier.
3. It's better to define vars and eqns as functions of n.

The following is my solution:

nmax = 5; T = 70; repeat = 1000;
r[n_] := RandomInteger[{1, n}];
vars[n_] := Table[x[j], {j, n}]; 
eqns[n_, ξ_] := 
 Table[x[j]'[t] == 
   x[j][t] (1 - x[j][t] - Sum[ξ[[j, k]] x[k][t] Boole[k != j], {k, n}]), {j, n}];

ans = Table[Module[{ics = {x[1][0] == 0.7}},
    Table[
     Module[{xi = RandomReal[1, {n, n}], sol, icvalue},
      sol = NDSolveValue[{eqns[n, xi], ics}, vars[n], {t, 0, T}];
      icvalue = Through@sol[T];
      ics = {Thread[Table[x[j][0], {j, n}] == icvalue], 
        x[n + 1][0] == 0.01 icvalue[[r@n]]};
      {sol, xi}], {n, nmax}]],
   {repeat}];