Thanks in advance for your help. I'm not looking for an explicit solution to my issue, but rather to be pointed in the right direction.
I have done all of the work to understand and recreate the model from scratch. I presented the qualitative part for a class project. What I am doing now is taking that project a step farther and trying to first solve the system in MATLAB with runge-kutta (or any method that works). Finally, I want to dive into the theory behind the numerical analysis to find out why the chosed method converges.
I have found that I can create a plot with roughly the same shape, but there are several problems:
1.) The time-scale over which the change occurs is three times that of the above plot.
2.) The range of function values is is vastly wrong.
3.) The desired shapes only occur if I tweak the initial conditions to be significantly different than what is shown near t=0 above.
So, finally, my question is simply to ask what might be causing the disparity I am seeing in my plots?
Code thus far:
r_b = 1.52;
k_b = 355;
alph = 1.11;
bet = 43200;
r_e = 0.92;
k_e = 1;
p = 0.00195;
r_s = 0.095;
k_s = 25440;
tspan = [0 200];
init = [1 1 1];
[t, Y] = ode45(@(t,y) odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s), tspan, init);
function dydt = odefcn(t, y, r_b, k_b, alph, bet, r_e, k_e, p, r_s, k_s)
dydt = [ r_b*y(1)*(1 - y(1)/(k_b*y(2))) - bet*(y(1)^2/((alph*y(2))^2 + y(1)^2));
r_s*y(2)*(1 - (y(2)*k_e)/(k_s*y(3)));
r_e*y(3)*(1 - (y(3)/k_e)) - p*y(1)/y(2)