Numerical Solutions For Unary Vaccine

Compound variable function Cp for carrying capacity
Several parameter sets at once

Scott’s Code

Pars is now a matrix of parameters. Each row is a single simulation run. I have set it up with three, but you can set it up to do however many you want. You just need to add rows to the Pars matrix and change the for loop to run through them all (the i<= in the loop). It will save all your plots in “Plots” and you can barf them out at will.

Pars={{(Subscript[b, v])->0.5,(Subscript[b, L])->300,(Subscript[k, v])->.0000000001,(Subscript[k, L])->.0000000001,(Subscript[r, u])->.034,(Subscript[r, v])->.027,τ->45, X->2800000000, Y -> 640000000},{(Subscript[b, v])->1,(Subscript[b, L])->300,(Subscript[k, v])->.0000000001,(Subscript[k, L])->.0000000001,(Subscript[r, u])->.034,(Subscript[r, v])->.027,τ->45, X->2800000000, Y -> 640000000},{(Subscript[b, v])->2,(Subscript[b, L])->300,(Subscript[k, v])->.0000000001,(Subscript[k, L])->.0000000001,(Subscript[r, u])->.034,(Subscript[r, v])->.027,τ->45, X->2800000000, Y -> 640000000}};

Module[{i},Plots={};For[i=1,i<=3,i++,Cu[t]=(  1- (((Subscript[H, u])[t]+  (Subscript[H, v])[t])/ X)   );
Cv[t]=(  1- (((Subscript[H, u])[t]+  (Subscript[H, v])[t])/ Y)   );Sols=NDSolve[

{
V'[t]==(-(Subscript[k, v])V[t] (Subscript[H, u])[t]+(Subscript[b, v])(Subscript[H, v])[t] /.Pars[[i]]), 

L'[t]==(((Subscript[b, L])(Subscript[k, L])L[t-τ](Subscript[H, u])[t-τ]-(Subscript[k, L])L[t](Subscript[H, u])[t])/.Pars[[i]]),

(Subscript[H, u])'[t]==((  (Subscript[r, u])(Subscript[H, u])[t]Cu[t]-(Subscript[k, v])V[t](Subscript[H, u])[t]  -(Subscript[k, L])L[t](Subscript[H, u])[t]       )/.Pars[[i]]),

(Subscript[H, v])'[t]==( ((Subscript[r, v])(Subscript[H, v])[t]Cv[t]+(Subscript[k, v])V[t](Subscript[H, u])[t])        /.Pars[[i]]),

(Subscript[H, v])[0]==0,L[t/;t≤ 0]==100000,V[0]==0,(Subscript[H, u])[t/;t≤0]==100000000,WhenEvent[t==80,V[t]-> 100000]},

{V[t],(Subscript[H, u])[t],(Subscript[H, v])[t],L[t]},{t,0,900}   ];Plots=Append[Plots,LogPlot[{Evaluate [(Subscript[H, u])[t]/. Sols],Evaluate [(Subscript[H, v])[t]/. Sols],Evaluate [V[t]/. Sols],Evaluate [L[t]/. Sols]},{t,0,400},PlotRange->{0.01,1000000000000},Frame->True,PlotStyle->{Red,Blue,Green,Gray},PlotLegends->{"(Subscript[H, u])","(Subscript[H, v])","V","L"}]]];Clear[t]];

Plots