read "C:/Users/nikit/.maplesoft/DMB.txt";
                          First Written: Nov. 2021 


This is DMB.txt, A Maple package to explore Dynamical models in Biology (both 

   discrete and continuous)


   accompanying the class Dynamical Models in Biology, Rutgers University. 

      Taught by Dr. Z. (Doron Zeilbeger) 



               The most current version is available on WWW at:

          http://sites.math.rutgers.edu/~zeilberg/tokhniot/DMB.txt .

           Please report all bugs to: DoronZeil at gmail dot com .


             For general help, and a list of the MAIN functions,

        type "Help();". For specific help type "Help(procedure_name);" 


                        ------------------------------

            For a list of the supporting functions type: Help1();

             For help with any of them type: Help(ProcedureName);


                        ------------------------------

  For a list of the functions that give examples of Discrete-time dynamical 

     systems (some famous), type: HelpDDM();


             For help with any of them type: Help(ProcedureName);


                        ------------------------------

 For a list of the functions continuous-time dynamical systems (some famous) 

    type: HelpCDM();


             For help with any of them type: Help(ProcedureName);


                        ------------------------------

with(DocumentTools:-Layout);
with(DocumentTools:-Components);
CutOff := proc(b, c, alpha, beta) local a, z; for a from 0.01 by 0.01 to 400 do
z := OrbF(AllenSIRg(a, b, c, alpha, beta, x, y), [x, y], [0.2, 0.8], 10000,
10001); if 1.0/10^5 < z[1][1] then RETURN(a); end if; end do; end proc;



IncAlphaBeta := proc(alpha1, alphainc, alpha2, beta1, betainc, beta2, b1, binc,
b2, c1, cinc, c2, xaxis) local i, G, DT, xml, alpha, b, c, beta, count, co, L; i
:= 0; count := 0; for alpha from alpha1 by alphainc to alpha2 do for beta from
beta1 by betainc to beta2 do for b from b1 by binc to b2 do for c from c1 by
cinc to c2 do count := count + 1; end do; end do; end do; end do; G :=
Matrix(count, 5); for alpha from alpha1 by alphainc to alpha2 do for beta from
beta1 by betainc to beta2 do for b from b1 by binc to b2 do for c from c1 by
cinc to c2 do i := i + 1; co := CutOff(b, c, alpha, beta); G[i, 1] := b; G[i, 2]
:= c; G[i, 3] := alpha; G[i, 4] := beta; G[i, 5] := co; end do; end do; end do;
end do; DT := DataTable(identity = "DataTable1", variable = 'G', columnheader,
columnnames = ["b", "c", "alpha", "beta", "a0"], columnwidths = [50, 50, 50,
100, 200, 200, 200], tooltip = "my data table"); xml :=
Worksheet(Group(Input(Textfield(DT)))); DocumentTools:-InsertContent(xml); if
xaxis = "b" then L := DeleteColumn(G, [2, 3, 4]); elif xaxis = "c" then L :=
DeleteColumn(G, [1, 3, 4]); elif xaxis = "alpha" then L := DeleteColumn(G, [1,
2, 4]); elif xaxis = "beta" then L := DeleteColumn(G, [1, 2, 3]); end if;
pointplot(L, labels = [xaxis, "a0"]); end proc;


IncABC := proc(INI, a1, a2, ah, b1, b2, bh, c1, c2, ch, F, PARS, VARS) local i,
G, r, res, DT, xml, pa, pb, pc, count, F1; i := 0; count := 0; for pa from a1 by
ah to a2 do for pb from b1 by bh to b2 do for pc from c1 by ch to c2 do count :=
count + 1; end do; end do; end do; G := Matrix(count, 7); for pa from a1 by ah
to a2 do for pb from b1 by bh to b2 do for pc from c1 by ch to c2 do i := i + 1;
r := evalf(pa/(pb + pc)); F1 := subs(PARS[1] = pa, F); F1 := subs(PARS[2] = pb,
F1); F1 := subs(PARS[3] = pc, F1); res := OrbF(F1, VARS, [0.2, 0.8], 10000,
10003); G[i, 1] := pa; G[i, 2] := pb; G[i, 3] := pc; G[i, 4] := r; G[i, 5] :=
res[1]; G[i, 6] := res[2]; G[i, 7] := res[3]; end do; end do; end do; DT :=
DataTable(identity = "DataTable1", variable = 'G', columnheader, columnnames =
["a", "b", "c", "R0", "n = 10000", "n = 10001", "n = 10002"], columnwidths =
[50, 50, 50, 100, 200, 200, 200], tooltip = "my data table"); xml :=
Worksheet(Group(Input(Textfield(DT)))); DocumentTools:-InsertContent(xml); end
proc;

eq3Allen := proc(a, b, c, x, y) local pt; pt := [y*(1 - exp(-a*x))/(b + c), 1 -
x*(1 + c/b)]; end proc;
fixBC := proc(a1, a2, inc, b, c, INI, alpha, beta) local i, L, F, f, a, holda,
R, M, T, xml, s; if 1 < b + c or not type(INI, list) then print('bad*input');
RETURN(FAIL); end if; s := 0; for i from a1 by inc to a2 do s := s + 1; end do;
M := Matrix(s, 4); i := 1; for holda from a1 by inc to a2 do F :=
AllenSIRg(holda, b, c, alpha, beta, x, y); L := OrbF(F, [x, y], INI, 1000,
1010); R := evalf(holda/(b + c)); M[i, 1] := holda; M[i, 2] := R; M[i, 3] :=
L[1]; f := eq3Allen(holda, b, c, L[-1][1], L[-1][2]); M[i, 4] := f; i := i + 1;
end do; T := DataTable(identity = "Table A", variable = 'M', columnheader,
columnnames = ["a", "R", "Numeric Equilibrium Solution at n = 1010", "Solution
According to eq3"], columnwidths = [50, 50, 300, 300], tooltip = "my data
table"); xml := Worksheet(Group(Input(Textfield(T))));
DocumentTools:-InsertContent(xml); end proc;
fixAC := proc(b1, b2, inc, a, c, INI, alpha, beta) local i, L, holdL, F, f,
holdf, b, holdb, R, holdR, M, T, xml, s; if 1 < b1 + c or 1 < b2 + c or not
type(INI, list) then print('bad*input'); RETURN(FAIL); end if; s := 0; for i
from b1 by inc to b2 do s := s + 1; end do; M := Matrix(s, 4); i := 1; for holdb
from b1 by inc to b2 do F := AllenSIRg(a, holdb, c, alpha, beta, x, y); L :=
OrbF(F, [x, y], INI, 1000, 1010); R := evalf(a/(holdb + c)); M[i, 1] := holdb;
M[i, 2] := R; M[i, 3] := L[1]; f := eq3Allen(a, holdb, c, L[-1][1], L[-1][2]);
M[i, 4] := f; i := i + 1; end do; T := DataTable(identity = "Table B", variable
= 'M', columnheader, columnnames = ["b", "R", "Numeric Equilibrium Solution at n
= 1010", "Solution According to eq3"], columnwidths = [50, 50, 300, 300],
tooltip = "my data table"); xml := Worksheet(Group(Input(Textfield(T))));
DocumentTools:-InsertContent(xml); end proc;
fixAB := proc(c1, c2, inc, a, b, INI, alpha, beta) local i, L, holdL, F, f,
holdf, c, holdc, R, holdR, s, T, xml, M; if 1 < c1 + b or 1 < c2 + b or not
type(INI, list) then print('bad*input'); RETURN(FAIL); end if; s := 0; for i
from c1 by inc to c2 do s := s + 1; end do; i := 1; M := Matrix(s, 4); for holdc
from c1 by inc to c2 do F := AllenSIRg(a, b, holdc, alpha, beta, x, y); L :=
OrbF(F, [x, y], INI, 1000, 1010); R := evalf(a/(holdc + b)); M[i, 1] := holdc;
M[i, 2] := R; M[i, 3] := L[-1]; f := eq3Allen(a, b, holdc, L[-1][1], L[-1][2]);
M[i, 4] := f; i := i + 1; end do; T := DataTable(identity = "Table C", variable
= 'M', columnheader, columnnames = ["c", "R", "Numeric Equilibrium Solution at n
= 1010", "Solution According to eq3"], columnwidths = [50, 50, 300, 300],
tooltip = "my data table"); xml := Worksheet(Group(Input(Textfield(T))));
DocumentTools:-InsertContent(xml); end proc;

AllenSI := proc(a, b, x) ToSys(2, x, x[1]*(1 - b) + (1 - x[1])*(1 -
exp(-a*x[2])))[1]; end proc;
AllenSI(a, b, x);
             [x[1] (1 - b) + (1 - x[1]) (1 - exp(-a x[2])), x[1]]

AllenSIVar := proc(a, b, x) ToSys(2, x, x[1]*(1 - b) + (1 - x[1])*(1 -
exp(-a*x[2])))[2]; end proc;
with(LinearAlgebra);

ConjTwo := proc(INI, a1, a2, ah, b1, b2, bh, F, PARS, VARS) local i, G, r, res,
DT, xml, pa, pb, count, F1; i := 0; count := 0; for pa from a1 by ah to a2 do
for pb from b1 by bh to b2 do count := count + 1; end do; end do; G :=
Matrix(count, 8); for pa from a1 by ah to a2 do for pb from b1 by bh to b2 do i
:= i + 1; r := evalf(pa/pb); F1 := subs(PARS[1] = pa, F); F1 := subs(PARS[2] =
pb, F1); F1 := subs(PARS[3] = 0, F1); res := OrbF(F1, VARS, INI, 1000, 1003);
G[i, 1] := pa; G[i, 2] := pb; G[i, 3] := 0; G[i, 4] := r; G[i, 5] := res[1, 1];
G[i, 6] := res[2, 1]; G[i, 7] := res[3, 1]; G[i, 8] := evalf(subs(x = res[1][1],
(1 - x)*(1 - exp(-pa*x))/pb)); end do; end do; DT := DataTable(identity =
"DataTable1", variable = 'G', columnheader, columnnames = ["a", "b", "c", "R0",
"n=1000", "n=1001", "n=1002", "SI Model Equation"], columnwidths = [50, 50, 50,
100, 150, 150, 150, 150], tooltip = "my data table"); xml :=
Worksheet(Group(Input(Textfield(DT)))); DocumentTools:-InsertContent(xml); end
proc;