Numerical solution is off for system of ODEs

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 been plugging away at solving a system of non-linear, first order ODEs in MATLAB. The system was solved numerically in this study: http://web.math.ku.dk/~moller/e04/bio/ludwig78.pdf

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:

% System Parameters:
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);
subplot(3,1,1);
plot(t,Y(:,1),'b');
title('Budworm Density');

subplot(3,1,2)
plot(t,Y(:,2),'g');
title('Branch Density');

subplot(3,1,3);
plot(t,Y(:,3),'r');
title('Foliage Condition');

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)
        ];

end

Answer 1

[SandeepKumar R]

To be consistent use years also in your ode function (better to stick to year). Secondly play with tolerances of ode( Reltol,AbsTol). If that doesn't work change the ode solver (23,15s..etc..)