#M11.txt: Maple code for Lecture 11 of Dynamical Models in Biology taught by Dr. Z.
Help11:=proc(): print(` SFPe(f,x), Orbk(k,z,f,INI,K1,K2) `):end:

#SFPe(f,x): The set of fixed points of x->f(x) done exactly (and allowing symbolic parameters), followed by the condition of stability (if it is netween -1 and 1 it is stable)
#Try: FPe(k*x*(1-x),x);
#VERSION OF Oct. 12, 2021 (avoiding division by 0)
SFPe:=proc(f,x) local f1,L,i,M:
f1:=normal(diff(f,x)):
L:=[solve(numer(f-x),x)]: 
M:=[]:

for i from 1 to nops(L) do
 if subs(x=L[i],denom(f1) )<>0 then
  M:=[op(M),[L[i],normal(subs(x=L[i],f1))]]:
 fi:
od:
M:

end:


#Added after class
#Orbk(k,z,f,INI,K1,K2): Given a positive integer k, a letter (symbol), z, an expression f of z[1], ..., z[k] (representing a multi-variable function of the variables z[1],...,z[k]
#a vector INI representing the initial values [x[1],..., x[k]], and (in applications) positive integres K1 and K2, outputs the
#values of the sequence starting at n=K1 and ending at n=K2. of the sequence satisfying the difference equation
##x[n]=f(x[n-1],x[n-2],..., x[n-k+1]):
#This is a generalization to  higher-order difference equation of procedure Orb(f,x,x0,K1,K2). For example
#Orbk(1,z,5/2*z[1]*(1-z[1]),[0.5],1000,1010);  should be the same as
#Orb(5/2*z[1]*(1-z[1]),z[1],[0,5],1000,1010);
#Try:
#Orbk(2,z,(5/4)*z[1]-(3/8)*z[2],[1,2],1000,1010);
Orbk:=proc(k,z,f,INI,K1,K2) local L, i,newguy:
L:=INI:  #We start out with the list of initial values

if not (type(k, integer) and type(z,symbol) and type(INI,list) and nops(INI)=k and type(K1,integer) and type(K2,integer) and K1>0 and K2>K1) then #checking that the input is OK
 print(`bad input`):
 RETURN(FAIL):
fi:

while nops(L)<K2 do
 newguy:=subs({seq(z[i]=L[-i],i=1..k)},f):   #Using what we know about the value yesterday, the day before yesterday, ... up to k days before yesterday we find the value of the sequence today
 L:=[op(L),newguy] :  #we append the new value to the running list of values of our sequence
od:

[op(K1..K2,L)]:

end:


####STAFT FROM M9.txt
#M9.txt: Maple Code for "Dynamical models in Biology" (Math 336) taught by Dr. Z., Lecture 9

Help9:=proc(): print(` Orb(f,x,x0,K1,K2), Orb2D(f,x,x0,K) , FP(f,x) , SFP(f,x) , Comp(f,x)  `):end:

#Orb(f,x,x0,K1,K2): Inputs an expression f in x (desccribing) a function of x, an initial point, x0, and a positive integer K, outputs 
#the values of x[n] from n=K1 to n=K2. Try: where x[n]=f(x[n-1]), . Try:
#Orb(2*x*(1-x),x,0.4,1000,2000);
Orb:=proc(f,x,x0,K1,K2) local x1,i,L:
x1:=x0:
for i from 1 to K1 do
 x1:=subs(x=x1,f):   #we don't record the first values of K1, since we are interested in the long-time behavior of the orbit
od:

L:=[x1]:

for i from K1 to K2 do
 x1:=subs(x=x1,f):    #we compute the next member of the orbit
 L:=[op(L),x1]: #we append it to the list
od:

L: #that's the output

end:

#Orb2D(f,x,x0,K): 2D version of Orb(f,x,x0,0,K), just for illustration
Orb2D:=proc(f,x,x0,K) local L,L1,i:
L:=Orb(f,x,x0,0,K):
L1:=[[L[1],0],[L[1],L[2]], [L[2],L[2]] ] :
for i from 3 to nops(L) do
 L1:=[op(L1),[L[i-1],L[i]],[L[i],L[i]]]:
od:
L1:
end:


#FP(f,x): The list of fixed points of the map x->f where f is an expression in x. Try:
#FP(2*x*(1-x),x);
FP:=proc(f,x)
evalf([solve(f=x)]):
end:


#SFP(f,x): The list of stable fixed points of the map x->f where f is an expression in x. Try:
#SFP(2*x*(1-x),x);
SFP:=proc(f,x) local L,i,f1,pt,Ls:
L:=FP(f,x):  #The list of fixed points (including complex ones)

Ls:=[]:         #Ls is the list of stable fixed points, that starts out as the empty list

f1:=diff(f,x):   #The derivative of the function f w.r.t. x

for i from 1 to nops(L) do
pt:=L[i]:

if abs(subs(x=pt,f1))<1 then

 Ls:=[op(Ls),pt]:   # if pt,  is stable we add it to the list of stable points

fi:

od:

Ls:    #The last line is the output

end:

#Comp(f,x): f(f(x))
Comp:=proc(f,x): normal(subs(x=f,f)):end: