1 years, 47 days

## System solutions repeating itself...

 >
 >
 (1)
 >
 (2)
 >
 >

@Carl Love thanks very much. I corrected that, but i still have the system solutions repeating itself. This #is what i got.

## This is what i observe and still need so...

@tomleslie thank you very much for the job well done. This is the problem i encounter
(1) My maple version is not running the graphs as it ought to be, but it ran through in yours

Maple Worksheet - Error

Failed to load the worksheet /maplenet/convert/someODEs2.mw .

Maple Worksheet - Error

Failed to load the worksheet /maplenet/convert/someODEs2.mw .

Maple Worksheet - Error

Failed to load the worksheet /maplenet/convert/someODEs2.mw .

Maple Worksheet - Error

Maple Worksheet - Error

.

. you will see it in the code i uploaded. i mean it was giving me something like
PLOT( )
PLOT( ) and so on.

(2) i want to obtain the numerical values  from the series code you wrote as

(3)Also, if the answers are thesame i think the combined graphs should fit together
i mean the code you wrote as
display([v1, a1]); display([v2, a2]); display([v4, a4])

## @Christian Wolinski  I initialized ...

@Christian Wolinski
I initialized at B(0) := .50; C(0) := .30; DD(0) := .21; E(0) := .14; F(0) := .70; G(0) := .45; H(0) := .14 but it wasnt running smoothly

## The first two points raised by you are v...

@mmcdara thanks for your help. The first two points raised by you are very correct. B(T) is the priori model of the empirical observations iV. So I want iV to be considered as a observation for C(T) or DD(T). The model is resumed by 16 empirical observations by 7 matrix dimensions. Also, I want to know if transmissibility parameters like beta1 and betao values can be fitted to the experimental observations iV. You can please go forward with the coding/computations. Thanks.

## There is no relationship between i(T) a...

@mmcdara
There is no relationshipt between i(t) and B(T), C(T),....H(T)....it is wrong. i intend to write either B(T) or....H(T) where i wrongly wrote i(T).

My intention is to fit the data values with each state variables and see the behavior on the two on the same graph

## Fitting real life assumed data and compa...

@mmcdara 777 thanks a lot for your expalantions and the code correction.
These are my intentions;
(1) The system of equations describes an epidemic transmission which i got an assumed data values for infected individuals per year in different localities earlier listed in the begining of the code. i intend to fit this data to the model equations which still retains its parameters and variable values.
I want to fit this data values with each variables to see the relationship/behavior between real life data and parameters/variable values involved in the model system equations.
Or if you have another fitting method you can help with that and i will gladly apprreciate.
Please i need help on that.
Thanks

## power series code not working for my ver...

@tomleslie
The code to obtain the power series solution to this system is not working, Please i need help on a new code that will support my maple version

 >
 (1)
 >
 >
 >
 >
 (2)
 >
 >

## complex or singularity error...

@tomleslie good day, thanks for your attention. What can cause complex or singularity error in this case as shown below

 >
 >
 >
 >
 >

## Variation of parameters...

@tomleslie Thanks once again for your good work and special attention to this program.
My question is
(1) You have varied the parameter "xi" at three different values. I checked the code for alteration so that i can also vary other parameters like "Mh, psi, mu" etc. but i couldnt get it.
How do i vary other parameters.

## Salient points raised on the maple codes...

@tomleslie thank you very much for taking your time to attend to my problems, i'm very grateful. Here are some points i need to raise;
(1)A particular code you wrote as
interface(rtablesize = [22, 10]); M := Matrix([`~`[lhs](ans(0)), seq(`~`[rhs](ans(j)), j = 0 .. 20)])
is not running on my maple sheet. I dont know if my version cant support that or better still, another code that can give the matrix is appreciated.
(2) The code i wrote for the series is the differential transformation method(DTM) series applied to model equations. My intentions were to see wether the series solution for DTM will agree with the numerical values/solutions using RK4
(3) Also, i need codes for parametric plots, e.g., for Mh, beta1 etc.
Once again, I appreciate @tomleslie

 >
 >
 >
 >
 >
 (1)
 > # # Generate power series solutions of order 12 # (Maple default is 6) for all ODEs #   Order:=12:   ans2:= convert~          ( evalf~            ( dsolve              ( {ODEs, bcs},                indets([ODEs], function(name)),                'series'              )            ),            polynom          ); # # Plot the "numeric" solution obtained earlier # and the equivalent power series solution # obtaine above on th same graph. # # Note that the range of the plots is reduced, # because even with 12-th order polynomials, # the power series soluton and the numerical # solution start to deviate for T>~6 #   for j in indets([ODEs], function(name)) do       display       ( [ odeplot           ( ans,             [T, j],             T = 0 .. 10,             thickness = 4, style = point,             color = red           ),           plot           ( eval(j, ans2),             T = 0 .. 10,             thickness = 2, style = point,             color=blue           )         ],         title = cat(convert(j,string), " (Numerical and PowerSeries compared)"),         titlefont = [tims, bold, 20]       )    end do
 >

 >
 >
 >
 >
 >
 (1)
 > # # Generate power series solutions of order 12 # (Maple default is 6) for all ODEs #   Order:=12:   ans2:= convert~          ( evalf~            ( dsolve              ( {ODEs, bcs},                indets([ODEs], function(name)),                'series'              )            ),            polynom          ); # # Plot the "numeric" solution obtained earlier # and the equivalent power series solution # obtaine above on th same graph. # # Note that the range of the plots is reduced, # because even with 12-th order polynomials, # the power series soluton and the numerical # solution start to deviate for T>~6 #   for j in indets([ODEs], function(name)) do       display       ( [ odeplot           ( ans,             [T, j],             T = 0 .. 10,             thickness = 4, style = point,             color = red           ),           plot           ( eval(j, ans2),             T = 0 .. 10,             thickness = 2, style = point,             color=blue           )         ],         title = cat(convert(j,string), " (Numerical and PowerSeries compared)"),         titlefont = [tims, bold, 20]       )    end do
 >

 >
 >
 >
 >
 >
 (1)
 > # # Generate power series solutions of order 12 # (Maple default is 6) for all ODEs #   Order:=12:   ans2:= convert~          ( evalf~            ( dsolve              ( {ODEs, bcs},                indets([ODEs], function(name)),                'series'              )            ),            polynom          ); # # Plot the "numeric" solution obtained earlier # and the equivalent power series solution # obtaine above on th same graph. # # Note that the range of the plots is reduced, # because even with 12-th order polynomials, # the power series soluton and the numerical # solution start to deviate for T>~6 #   for j in indets([ODEs], function(name)) do       display       ( [ odeplot           ( ans,             [T, j],             T = 0 .. 10,             thickness = 4, style = point,             color = red           ),           plot           ( eval(j, ans2),             T = 0 .. 10,             thickness = 2, style = point,             color=blue           )         ],         title = cat(convert(j,string), " (Numerical and PowerSeries compared)"),         titlefont = [tims, bold, 20]       )    end do
 >

 >
 >
 >
 >
 >
 (1)
 > # # Generate power series solutions of order 12 # (Maple default is 6) for all ODEs #   Order:=12:   ans2:= convert~          ( evalf~            ( dsolve              ( {ODEs, bcs},                indets([ODEs], function(name)),                'series'              )            ),            polynom          ); # # Plot the "numeric" solution obtained earlier # and the equivalent power series solution # obtaine above on th same graph. # # Note that the range of the plots is reduced, # because even with 12-th order polynomials, # the power series soluton and the numerical # solution start to deviate for T>~6 #   for j in indets([ODEs], function(name)) do       display       ( [ odeplot           ( ans,             [T, j],             T = 0 .. 10,             thickness = 4, style = point,             color = red           ),           plot           ( eval(j, ans2),             T = 0 .. 10,             thickness = 2, style = point,             color=blue           )         ],         title = cat(convert(j,string), " (Numerical and PowerSeries compared)"),         titlefont = [tims, bold, 20]       )    end do
 >

 >
 >
 >
 >
 >
 (1)
 > # # Generate power series solutions of order 12 # (Maple default is 6) for all ODEs #   Order:=12:   ans2:= convert~          ( evalf~            ( dsolve              ( {ODEs, bcs},                indets([ODEs], function(name)),                'series'              )            ),            polynom          ); # # Plot the "numeric" solution obtained earlier # and the equivalent power series solution # obtaine above on th same graph. # # Note that the range of the plots is reduced, # because even with 12-th order polynomials, # the power series soluton and the numerical # solution start to deviate for T>~6 #   for j in indets([ODEs], function(name)) do       display       ( [ odeplot           ( ans,             [T, j],             T = 0 .. 10,             thickness = 4, style = point,             color = red           ),           plot           ( eval(j, ans2),             T = 0 .. 10,             thickness = 2, style = point,             color=blue           )         ],         title = cat(convert(j,string), " (Numerical and PowerSeries compared)"),         titlefont = [tims, bold, 20]       )    end do
 >

 >
 >
 >
 >
 >
 (1)
 > # # Generate power series solutions of order 12 # (Maple default is 6) for all ODEs #   Order:=12:   ans2:= convert~          ( evalf~            ( dsolve              ( {ODEs, bcs},                indets([ODEs], function(name)),                'series'              )            ),            polynom          ); # # Plot the "numeric" solution obtained earlier # and the equivalent power series solution # obtaine above on th same graph. # # Note that the range of the plots is reduced, # because even with 12-th order polynomials, # the power series soluton and the numerical # solution start to deviate for T>~6 #   for j in indets([ODEs], function(name)) do       display       ( [ odeplot           ( ans,             [T, j],             T = 0 .. 10,             thickness = 4, style = point,             color = red           ),           plot           ( eval(j, ans2),             T = 0 .. 10,             thickness = 2, style = point,             color=blue           )         ],         title = cat(convert(j,string), " (Numerical and PowerSeries compared)"),         titlefont = [tims, bold, 20]       )    end do
 >

 >
 >
 >

## to modify the previous code to mine...

 > restart: # # Statement not terminated in the following - fixed it #   B[0]:= 1.7e8:  C[0]:=0:  P[0]:=0.1: E[0]:=2: # # Why does this (commented) Digits command exist? # #Digits:=10: # # OP never uses 'gamma' in this worksheet - why the # hell unprotect an inbuilt variable for no purpose??? # unprotect('gamma'):  A__h := 10:  beta__1 := 0.34e-1:  psi__o := 0.25e-1:  mu := 0.4e-3:  theta := .7902:  alpha := .11:  k := 0.136e-3:  z := 0.5e-1:  delta := .7:  eta := .134:  xi := .12:  beta := 0.1e-1:  Q1 := 1:  Q2 := 1:  h := .1: tau:=1.0;   kappa:=2.0; # # The variables 'tau' and 'kappa' are used in the following, # but neither of these is defined. Gave them a "random" # values here just to assist debug #
 (1)
 > #Step 1 n:=100: # # Fixed the names of all of the following for some kind # of consistency # lambda1[n] := 0: lambda2[n] := 0: lambda3[n] := 0: lambda4[n] := 0: u0[0] := 0:
 > for i from 0 to n-1 do      B[i+1] := h*(theta*E[i]+A__h+B[i])/(1+ psi__o*h*P[i])/(P[i]+xi)+mu+u0[i]+beta*C[i]; # # Following statement is a recursive assignment. # It was a recursive assignment the last time the # OP posted some code, and it is still a recursive # assignment now!!!! # # Let's see if I can make this explanation simple enough # for the brain-dead. # # C[i+1] is undefined. Therefore, assigning C[i+1] in # terms of C[i+1] is doomed to failure. # # Changing C[i+1] to C[i] on the rhs is *may* be incorrect # but at least the statement can now be evaluated!! # #     C[i+1] := h*B[i+1](beta*C[i]-psi[o]*P[i])/(P[i]+xi)-(alpha+mu+k)*C[i+1]; # # I'm also assumin that psi[o] ought to be psi__o, so # I made this change #      C[i+1] := h*B[i+1](beta*C[i]-psi__o*P[i])/(P[i]+xi)-(alpha+mu+k)*C[i];      P[i+1] := h*z*C[i+1]/(1+delta+eta);      E[i+1] := h*(alpha*C[i+1]+B[i+1]*u0[i]+E[i])/(1+theta+mu); # # The following statement has one more oopening parenthesis (ie '(') # than closing parenthesis - again basic syntax # # I randomly added a closing parenthesis, just to make it syntactically # correct, but the odds that I got this closing parenthesis in the # desired location are almost nil!!!!! # #   lambda1[n-i-1] := -h*(Q1+(beta*C[i]+(psi__o*P[i]/(P[i]+xi)))*lambda2[n-i]-lambda1[n-i]-mu+u0[i]/(1+beta*C        [i]+P[i]/(P[i]+xi)); #   lambda1[n-i-1] := -h*(Q1+(beta*C[i]+(psi__o*P[i]/(P[i]+xi)))*lambda2[n-i]-lambda1[n-i]-mu+u0[i])/(1+beta*C[i]+P[i]/(P[i]+xi)); # # The variable 'kappa' in the following is undefined #   lambda2[n-i-1] := -h*(-beta*B[i+1]*lambda1[n-i-1]+Q2-alpha-mu-kappa+z*lambda3[n-i]+alpha*lambda4[n-i]-lambda2[n-i])/(-beta*B[i+1]+1);   lambda3[n-i-1] := h*((psi__o*B[i+1]*P[i+1]/(P[i+1]+xi)^2+theta)*lambda1[n-i-1]-lambda2[n-i-1]*psi__o*B[i+1]*P[i+1]/(P[i+1]+xi)^2-lambda3[n-i])/(1+delta+eta);   lambda4[n-i-1] := h*lambda4[n-i]/(1-theta-mu); # # The variable 'tau' was undefined, so (obviously) the # following cannot be fully evaluated #      R1:= (lambda1[n-i-1]-lambda4[n-i-1])*B[i+1]/tau; # # Since R1 has not been fully evaluated, u0[i+1] cannot # be fully evaluated either - and so the chaos continues! #            u0[i+1] := min(1, max(R1, 0));         end do:
 > # # These two command are assigned to nothing. They are # are never used for anything - why do they exist?? # seq(i,i=0..30): seq(B[i],i=0..30):
 > with(plots): with(DEtools): # # The quantity 'h', was earlier defined to be 0.1 # Thus h(t)=0.1, h(0)=0.1 and so on. Probably not # what the OP desired. Even funnier is when the # OP uses initial conditions such as h(0)=400, # which of course will evaluate to 0.1=400. This # is a whole new level of misunderstanding # # Obviously this mean that the definition of these # ODEs is completely incorrect - so it is a complete\ # waste of time trying to solve them! # DE1 := diff(a(t),t) = A__h-beta__1*a(t)*g(t)-psi__o*a(t)*v(t)/(v(t)+xi)-mu*a(t)+theta*m(t); DE2 := diff(g(t),t) = beta__1*a(t)*g(t)-psi__o*a(t)*v(t)/(v(t)+xi)-(alpha+mu+k)*g(t); DE3 := diff(v(t),t) = z*g(t)-(delta+eta)*v(t); DE4 := diff(m(t),t) = alpha*g(t)-(theta+mu)*m(t); SIRsys := [DE1, DE2, DE3, DE4]; # # The commands DE1plot(), DE2plot(), DE3plot() DE4plot() # do not exist in Maple so the following will never # evaluate. OP (probably??) means DEplot(). However with # the incorrect use of the name 'h' described above, # even a syntactically correct command will not # produce anything meaningful. Correct command with # crap arguments will still fail!! # aplot := DEplot(SIRsys, [a(t), g(t), v(t), m(t)], t = 0 .. 100, a = 0.167e9 .. 0.17e9, [[a(0) = 0.17e9, g(0) = 0, v(0) = 4, m(0) = 2]], scene = [t, a(t)], thickness = 2, linecolor = red, stepsize = .1); gplot := DEplot(SIRsys, [a(t), g(t), v(t), m(t)], t = 0 .. 100, g = 0 .. 0.211, [[a(0) = 0.17e9, g(0) = 0, v(0) = 4, m(0) = 2]], scene = [t, g(t)], thickness = 2, linecolor = red, stepsize = .1): vplot := DEplot(SIRsys, [a(t), g(t), v(t), m(t)], t = 0 .. 100, v = 0 .. 0.1e4, [[a(0) = 0.17e9, g(0) = 0, v(0) = 4, m(0) = 2]], scene = [t, v(t)], thickness = 2, linecolor = red, stepsize = .1): mplot := DEplot(SIRsys, [a(t), g(t), v(t), m(t)], t = 0 .. 100, m = 0 .. 0.1e4, [[a(0) = 0.17e9, g(0) = 0, v(0) = 4, m(0) = 2]], scene = [t, m(t)], thickness = 2, linecolor = red, stepsize = .1): # # Did no syntax checking beyond this point becuase I got # really, REALLY BORED # L:=<|>: #points:= {seq(B[i],i=0..100)};
 Warning, plot may be incomplete, the following errors(s) were issued:    cannot evaluate the solution further right of .33370860e-5, probably a singularity
 Warning, plot may be incomplete, the following errors(s) were issued:    cannot evaluate the solution further right of .33370860e-5, probably a singularity Warning, plot may be incomplete, the following errors(s) were issued:    cannot evaluate the solution further right of .33370860e-5, probably a singularity Warning, plot may be incomplete, the following errors(s) were issued:    cannot evaluate the solution further right of .33370860e-5, probably a singularity
 > SIRsys; a(0);g(0);v(0); m(0); dsolve([SIRsys[], a(0) = 0.17e9, g(0) = 0, v(0) = 4, m(0) = 2]);
 (2)
 > p1:=pointplot(L, color=blue, symbol=diamond): display([p1,aplot]);
 >
 >

## I am trying to compare a code with my co...

@tomleslie thanks very much for being patient with me( a maple beginner). i am sending a code that ran successfully please compare with mine for futher assistance

 >
 > restart: H[0]:= 1.7e8:  J[0]:=0:  V[0]:=400: #Digits:=10: unprotect('gamma'); lambda:=5e5: mu:=0.003: beta:=4e-10: delta:=0: alpha:=0.043: sigma:=alpha+delta: k:=6.24: gamma:=0.65: A1:=1: A2:=1: h:=1:
 > #Step 1 n:=100:  lambda__1[n]:= 0:   lambda__2[n]:= 0:   lambda__3[n]:= 0: u__1[0]:= 0:  u__2[0]:= 0:
 > for i from 0 to n-1 do      H[i+1]:=(H[i]+delta*J[i]*h+lambda*h)/(1+mu*h+beta*V[i]*h*(1-u__1[i]));      J[i+1]:= (J[i]+beta*H[i+1]*V[i]*h*(1-u__1[i]))/(sigma*h+1);      V[i+1]:=(V[i]+k*J[i+1]*h*(1-u__2[i]))/(1+gamma*h); lambda__1[n-i-1]:= -(-lambda__1[n-i]-lambda__2[n-i]*beta*V[i+1]*h+lambda__2[n-i]*beta*V[i+1]*h*u__1[i]-h)/(mu*h+beta*V[i+1]*h-beta*V[i+1]*h*u__1[i]+1);      lambda__2[n-i-1]:= -(-lambda__2[n-i]-lambda__1[n-i-1]*delta*h-lambda__3[n-i]*k*h+lambda__3[n-i]*k*h*u__2[i])/(1+sigma*h);      lambda__3[n-i-1]:= (lambda__3[n-i]-lambda__1[n-i-1]*beta*H[i+1]*h+lambda__1[n-i-1]*beta*H[i+1]*h*u__1[i]+lambda__2[n-i-1]*beta*H[i+1]*h-lambda__2[n-i-1]*beta*H[i+1]*h*u__1[i])/(1+gamma*h);      R1:= (1/A1)*(lambda__1[n-i-1]-lambda__2[n-i-1])*beta*V[i+1]*H[i+1];      R2:= -(1/A2)*lambda__3[n-i-1]*k*J[i+1];      u__1[i+1]:= min(1, max(R1,0));      u__2[i+1]:= min(1, max(R2,0));     end do:
 > seq(i,i=0..30): seq(H[i],i=0..30):
 > with(plots): with(DEtools): DE1 := diff(s(t),t)=lambda-mu*s(t)-beta*s(t)*r(t)+delta*g(t); DE2 := diff(g(t),t)=beta*s(t)*r(t)-sigma*g(t); DE3 := diff(r(t),t)=k*g(t)-gamma*r(t); SIRsys:= [DE1, DE2, DE3]; splot:=DEplot(SIRsys, [s(t),g(t),r(t)], t = 0..100,s=1.67e8..1.7e8,[[s(0)=1.7e8,g(0)=0,r(0)=400]], scene=[t,s(t)], thickness=2, linecolor=red, stepsize=0.1): gplot:=DEplot(SIRsys, [s(t),g(t),r(t)], t = 0..100,g=0..400, [[s(0)=1.7e8,g(0)=0,r(0)=400]], scene=[t,g(t)], thickness=2, linecolor=red, stepsize=0.1): rplot:=DEplot(SIRsys, [s(t),g(t),r(t)], t = 0..100,r=0..1e3,[[s(0)=1.7e8,g(0)=0,r(0)=400]], scene=[t,r(t)], thickness=2, linecolor=red, stepsize=0.1): A:=<|>; #points:= {seq(H[i],i=0..100)};
 (1)
 > p1:=pointplot(A, color=blue, symbol=diamond): display([p1,splot]);
 > B:=<|>; #points:= {seq(J[i],i=0..100)}; p2:=pointplot(B, color=blue, symbol=soliddiamond): display([p2,gplot]);
 > C:=<|>; #points:= {seq(V[i],i=0..100)}; p3:=pointplot(C, color=blue, symbol=diamond): display([p3,rplot]);
 > DD:=<|>; #points:= {seq(u__1[i],i=0..100)}; pointplot(DD, color=black, symbol=diamond);
 > E:=<|>; #points:= {seq(u__2[i],i=0..100)}; pointplot(E, color=blue, symbol=diamond);
 >
 >
 >
 >
 >
 >
 >
 >
 >
 >
 (2)
 >

## I still encounter little errors constrai...

#Great job done by @tomleslie4564, thank you.

 > restart: B[0]:= 1.7e8:  C[0]:=0:  P[0]:=400: E[0]:=20 #Digits:=10: unprotect('gamma'):  A__h := 10:  beta__1 := 0.34e-1:  psi__o := 0.25e-1: mu := 0.4e-3: theta := .7902:  alpha := .11:  k := 0.136e-3:  z := 0.5e-1: delta := .7: eta := .134:  xi := .12: beta := 0.1e-1: Q1 := 1: Q2 := 1: h := .1:
 > #Step 1 n:=100: lambda__1[n] := 0: lambda__2[n] := 0: lambda__3[n] := 0: lambda__4[n] := 0: u__o[0] := 0:
 > for i from 0 to n-1 do      B[i+1] := h*(theta*E[i]+A__h+B[i])/(1+ psi__o*h*P[i])/(P[i]+xi)+mu+u__[o][i]+beta*C[i];      C[i+1] := h*B[i+1](beta*C[i]-psi[o]*P[i])/(P[i]+xi)-(alpha+mu+k)*C[i+1];      P[i+1] := h*z*C[i+1]/(1+delta+eta);      E[i+1] := h*(alpha*C[i+1]+B[i+1]*u__[o][i]+E[i])/(1+theta+mu);      lambda__1[n-i-1] := -h*(Q1+(beta*C[i]+(psi__o*P[i]/(P[i]+xi)))*lambda__2[n-i]-lambda__1[n-i]-mu+u__o[i]/(1+beta*C        [i]+P[i]/(P[i]+xi));      lambda__2[n-i-1] := -h*(-beta*B[i+1]*lambda__1[n-i-1]+Q2-alpha-mu-kappa+z*lambda__3[n-i]+alpha*lambda__4[n-i]-           lambda__2[n-i])/(-beta*B[i+1]+1);      lambda__3[n-i-1] := h*((psi__o*B[i+1]*P[i+1]/(P[i+1]+xi)^2+theta)*lambda__1[n-i-1]-lambda__2[n-i-1]*psi__o*B[i+1]*P      [i+1]/(P[i+1]+xi)^2-lambda__3[n-i])/(1+delta+eta);      lambda__4[n-i-1] := h*lambda__4[n-i]/(1-theta-mu);      R1:= (lambda__1[n-i-1]-lambda__4[n-i-1])*B[i+1]/tau;            u__o[i+1] := min(1, max(R1, 0));          end do:
 > seq(i,i=0..30): seq(B[i],i=0..30):
 > with(plots): with(DEtools): DE1 := diff(a(t),t) = A__h-beta__1*a(t)*g(t)-psi__o*a(t)*h(t)/(h(t)+xi)-mu*a(t)+theta*m(t); DE2 := diff(g(t),t) = beta__1*a(t)*g(t)-psi__o*a(t)*h(t)/(h(t)+xi)-(alpha+mu+k)*g(t); DE3 := diff(h(t),t) = z*g(t)-(delta+eta)*h(t); DE4 := diff(m(t),t) = alpha*g(t)-(theta+mu)*m(t); SIRsys := [DE1, DE2, DE3, DE4]; aplot := DE1plot(SIRsys, [a(t), g(t), h(t), m(t)], t = 0 .. 100, a = 0.167e9 .. 0.17e9, [[a(0) = 0.17e9, g(0) = 0, h(0) = 400, m(0) = 20]], scene = [t, a(t)], thickness = 2, linecolor = red, stepsize = .1): gplot := DE2plot(SIRsys, [a(t), g(t), h(t), m(t)], t = 0 .. 100, g = 0 .. 400, [[a(0) = 0.17e9, g(0) = 0, h(0) = 400, m(0) = 20]], scene = [t, g(t)], thickness = 2, linecolor = red, stepsize = .1): hplot := DE3plot(SIRsys, [a(t), g(t), h(t), m(t)], t = 0 .. 100, h = 0 .. 0.1e4, [[a(0) = 0.17e9, g(0) = 0, h(0) = 400, m(0) = 20]], scene = [t, h(t)], thickness = 2, linecolor = red, stepsize = .1): mplot := DE4plot(SIRsys, [a(t), g(t), h(t), m(t)], t = 0 .. 100, m = 0 .. 0.1e4, [[a(0) = 0.17e9, g(0) = 0, h(0) = 400, m(0) = 20]], scene = [t, m(t)], thickness = 2, linecolor = red, stepsize = .1): L:=<|>: #points:= {seq(B[i],i=0..100)};
 (1)
 > p1:=pointplot(L, color=blue, symbol=diamond): display([p1,aplot]);
 >
 >