#Dynamical Models in Bio Project 3

#Proj3.txt: Maple package for Project 3 of Math 336 (Rutgers Univ., Fall 2025 taught by Dr. Z. (aka Professor Doron Zeilberger))
read `DMB.txt`:


print(`This is Proj3.txt, a Maple package accompanying Project 3 of Dr. Z.'s Math 336 (Dynamical Models in Biology), Fall 2025`):
print(` Project Leader: Aurelia Altzman `):
print(`Other members: Adriana Ferreira, Palash Keswani`):
print(`For a list of the procedures done by this team type: Help3(); For Help with a specific procedure, type: Help3(ProcedureName);`):

with(plots):

Help3:=proc()
  if nargs=0 then
    print(`The available procedures are:  `):
    print(`F, f, OrbF, showExample, detectPeriod, showBiforcations`) 
  elif args[1] = detectPeriod then 
    print(`guesses the period of an orbit. Chaotic orbits should be at least 65 elements long. Returns 64 if the period is longer than 64 or does not exist`):
  elif args[1] = f    then print(`An example logistic equation. f(k) = k * x * (1-x)`):
  elif args[1] = OrbF then print(`Returns an Orbit of F(n) of length 17 w/ initial value 0.2`):
  elif args[1] = F    then print(`Logistic equation. Usage: F([2,3,4], k, x) = x - (x-2)(x-3)(x-4)/k`)
  elif args[1] = showExample then print(`Shows off the estimated and real stable values for varying values of k`):
                                  print(`Usage: showExample(F([3,4,5], k, x)`):
				  print(`(Output may be long!)`):
  elif args[1] = showBiforcations then print(`Shows a graph of the period doubling behavior of the logistic equation.`)  
  else
    print(`There is no Help for`, args):
  end:

end:
detectPeriod := proc(L) 
  if L[1] = L[2] then return 1 
  elif L[1] = L[3] then return 2
  elif L[1] = L[4] then return 3 
  elif L[1] = L[5] then return 4 
  elif L[1] = L[9] then return 8 
  elif L[1] = L[17] then return 16 
  elif L[1] = L[33] then return 32 
  elif L[1] = L[65] then return 64
  else return 64
  end:
end:
f := proc(k) return k*x*(1-x) end:
OrbF := proc(n) Orb(f(n), x, 0.2, 1000, 1064) end:
take := proc(L) L[1..detectPeriod(L)] end:
F:=proc(A,k,x) local i: x-mul(x-A[i],i=1..nops(A))/k:end: 
take := proc(L) L[1..detectPeriod(L)] end:
showBiforcations := proc() 
  Digits := 20:
  local L := []:
  local i:
  for i from 1 to 200 do 
    local n := 2 + 0.01 * i;
    local elem:
    for elem in take(OrbF(n)) do 
      L:= [op(L), [n, elem]]:
    end:
  end:
  pointplot(L);
  #return L 
end:

toSet:=proc(l) 
  {op(l)} 
end:
guessBiforcatedValues:=proc(f, x) 
	local s := subs(x=f, f): #; print(s) ;
	{op(SSS(s, x))}
end:
trueBiforcatedValues := proc(f, x) 
	local o := Orb(f, x, 0.5, 1000,1064):
	{op(take(o))} 
end:

example0 :=F([5,6,7], k, x) :
example2 := F([2,4,8],k, x):
#example1 := F([1,10,100], k, x):

showExample:=proc(f, k, x)
    local idx, s, guess;
    for idx from 0.1 by 0.01 to 100 do s := subs(k = idx, f); guess := guessBiforcatedValues(s, x); if guess <> {} then print(idx, guess); end if; if SSS(s, x) <> [] then print(SSS(s, x)); break; end if; end do;
end proc:
example3 := F([3,5,7,9], k, x) :
example4 := F([3, 8, 11], k, x):
example5 := F([1,2,7,8], k, x):
example6 := F([9,10,11], k, x):
example7 := F([8,9,10], k, x): 
example8 := F([8,9,10,11], k, x):
example9 := F([8,9,10,11,12],k,x):
example10 := F([2,3], k, x):