> |
restart:
with(plots):
with(plottools):
with(DEtools):
eqn1 := diff(V(t), t) = pi*p - (alpha + mu)*V(t),
V(0) = ic1:
eqn2 := diff(S(t), t) = alpha*V(t) + (1 - p)*pi - beta*S(t)*In(t)/N - mu*S(t),
S(0) = ic2:
eqn3 := diff(In(t), t) = beta*S(t)*In(t)/N - (mu + delta + localgamma)*In(t),
In(0) = ic3:
eqn4 := diff(R(t), t) = localgamma*In(t) - mu*R(t),
R(0) = ic4:
pi := 487845:
pVals := [0.2, 0.5, 0.7, 0.8]:
alpha := 0.054:
beta := 0.955:
mu := 0.005:
delta := 0.03:
localgamma := 0.935:
ic1 := 484465:
ic2 := 31999760:
ic3 := 26305:
ic4 := 12470:
N:=10000:
cols:=[red, blue, green, black]:
for k from 1 by 1 to numelems(pVals) do
p:=pVals[k]:
dsol := dsolve([eqn1, eqn2, eqn3, eqn4], numeric):
plt[k]:=odeplot( dsol,
[t, rhs(eqn1[1])],
t=0..200,
color=cols[k]
):
od:
display( [seq( plt[j], j=1.. numelems(pVals))],
legend=[seq( 'p'=pVals[j], j=1..numelems(pVals))],
labels=[ typeset(t), typeset(lhs(eqn1[1])) ],
labelfont=[times, bold, 16]
);
|