# Adding extra species in a model

I am currently working on a model of competition for resources between different species. I have 6 species competing for 3 resources. The first 3 species start at =0. the 4th species is added at t=1000 days, the 5th species is added at t=2000 days and the 6th species at t=5000. The first part is going well, but I can't seem to work out how I can get the last 3 species to be added later in the simulation. Can someone help me with this?

My code:

    (* Basic competition model *)
(* Ns[i] \[Rule] population abundance of species i *)
(* R[j] \[Rule] availability of resource j *)

(* constant equation stuff: *)
(* growth function of species i considering R of each resource at \
time t *)
\[Mu][i_, t_] :=
Min[(r[[#]]*R[#, t])/(Ks[[#, i]] + R[#, t]) & /@ Range[k]];

(* differential equations for Ns and R *)
dNs[i_] := D[Ns[i, t], t] == Ns[i, t]*(\[Mu][i, t] - m[[i]]);
dR[j_] :=
D[R[j, t], t] ==
d*(S[[j]] - R[j, t]) -
Sum[c[[j, i]]*\[Mu][i, t]*Ns[i, t], {i, 1, n}];
(* Initialize parameters: 5 species, 5 resources, chaos *)
n = 6; (* number of species *)
k = 3; (* number of resources *)
tMax = 10000; (* time steps (days) *)
r = ConstantArray[1, n]; (* maximum growth rate of species *)
m = ConstantArray[0.25, n]; (* mortality rate of species *)
S = {6, 10, 14}; (* supply concentration of resources *)
d = 0.25; (* system turnover rate *)
Ks = ( {
{1, 0.9, 0.3, 1.04, 0.34, 0.77},
{0.3, 1, 0.9, 0.71, 1.02, 0.76},
{0.9, 0.3, 1, 0.46, 0.34, 1.07}
});  (* Subscript[K, ji ]= half saturation constant of resource j \
for species i *)
c = ({
{0.04, 0.07, 0.04, 0.1, 0.03, 0.02},
{0.08, 0.08, 0.1, 0.1, 0.05, 0.17},
{0.14, 0.1, 0.1, 0.16, 0.06, 0.14}
});  (* Subscript[c, ji] = content of resource j in species i *)
(* Ns and Rs at time t = 0 *)
initNs = (Ns[#, 0] == 0.1 + #/100) & /@ Range[n];

initR = (R[#, 0] == S[[#]]) & /@ Range[k];

(* make a list of n many Ns equations, and k many R equations *)
sysNs = Array[dNs[#] &, n];
sysR = Array[dR[#] &, k];
(* make lists of variables to be solved for *)
resultsNs = Array[Ns[#, t] &, n];
resultsRs = Array[R[#, t] &, n];
sol2 = NDSolveValue[{sysNs, sysR, initNs, initR}, {resultsNs,
resultsRs}, {t, 0, tMax}];
solNs2 = sol2[[1]];
solR2 = sol2[[2]];
Plot[solNs2, {t, 0, 10000}]

• Would WhenEvent help? sol2 = NDSolveValue[{sysNs, sysR, initNs, initR, WhenEvent[t > 1000, Ns[4, t] -> 100]}, {resultsNs, resultsRs}, {t, 0, tMax}] – Ruud3.1415 Oct 2 at 10:46