This code runs except the graph and results are incredibly wrong

This code is running except the graph and the results are really wrong. I will provide the correct graph as a refference. I have checked all the equations and parameters but I cant seem to find why the answers are going to values that are *10^26

close all, clear, clc
load('EnvironmentalForcing.mat')
Bmax = 1;
uL_min = 1/6;
uI = 1/10;
e = 0.001;
Ap = 5000;
Pi = 930.27249;
Si = Pi/Ap;
Li = 0.01*Si;
Ii = 0;
Ri = uI*Ii;
Bi = 1e-4;
T_beta = zeros(1, length(tspan));
uL = zeros(1, length(tspan));
for i=1:length(tspan)
    if (T(i) <= 0) || (T(i)>=35)
        T_beta(i) = 0;
    else
        T_beta(i) = 0.00241*T(i)^2.06737 * (35-T(i))^0.72859;
    end 
end
for i=1:length(tspan)
    uL(i) = 1/sum(T_beta(1:i));
    j = 0;
    while 1/uL(i) > 1/uL_min
        j = j+1;
        uL(i) = 1/sum(T_beta(1+j: i));
    end
end
% Te = -0.35068 + 0.10789.*T -0.00214.*T.^2;
% for i = 1:length(T)
%     if T(i)>0 && T(i)<35
%         Tb(i) = (0.000241*(T(i)^2.06737))*((35-T(i))^0.72859);
%     end
% end
% j = 1;
% for i=1:length(t)
%     uL(i) = sum(Tb(j:i));
%     while uL(i) > uL_min
%         j = j+1;
%         uL(i) = sum(Tb(j:i));
%     end
% end
% B = Bmax*Tb;
% p = [Bmax, uL, uI, e, Te];
y0 = [Pi; Si; Li; Ii; Bi];
[t,y] = rk4(@PSLIR, tspan, y0, Bmax, uL, uI, e, T, tspan, Ap);
figure
hold on
plot(t,y(1,:),'k')
plot(t,y(2,:),'b')
plot(t,y(3,:),'g')
plot(t,y(4,:),'y')
plot(t,y(5,:),'m')
legend("S","L","I","R","P")

%%      Functions
function [dydt] = PSLIR(t_index, y0, Bmax, uL, uI, e, T, day, Ap)
if (isa(t_index, 'int8')) &&(t_index >= 2)
    t_index_rounded = round(t_index);
elseif (t_index <2)
    t_index_rounded = 1;
else
    t_index_rounded = round(t_index)-1;
end
if (t_index_rounded > 0) && (t_index_rounded <= length(T))
    Te = -0.35068 + 0.10789.*T(t_index_rounded) -0.00214.*T(t_index_rounded).^2;
    if (T(t_index_rounded) <=0) || (T(t_index_rounded) >=35)
        T_beta = 0;
    else
        T_beta = (0.000241*(T(t_index_rounded)^2.06737))*((35-T(t_index_rounded))^0.72859);
    end
    B = Bmax.*T_beta;
else
    error('Cant compute')
end
%assign variables
P = y0(1);
S = y0(2);
L = y0(3);
I = y0(4);
%R = y0(5);
Pb = y0(5);
dPbdt = (0.1724.*Pb-0.0000212.*Pb^2).*Te;
dPdt = ((1.33.*day(t_index_rounded)).*Te)+dPbdt;
dSdt = (-B.*S.*I)+dPdt./Ap;
dLdt = (B.*S.*I)-((uL(t_index_rounded).^-1).*L)+e;
dIdt = ((uL(t_index_rounded).^-1).*L)-((uI.^-1).*I);
dRdt = (uI.^-1).*I;
dydt = [dPdt; dSdt; dLdt; dIdt; dRdt];
end
function [t, y] = rk4(odeFunc, tspan, y0, Bmax, uL, uI, e, T, day, Ap)
%      N = length(tspan);
%     q = length(y0);
%     t0 = tspan(2);
%     h = tspan(2) - tspan(1);
%     t = zeros(N, 1);
%     y = zeros(q, N);
%     t(1) = t0;
%     y(:, 1) = y0;
N = length(tspan);
t = tspan;
y = zeros(length(y0),N);
y(:,1) = y0(:);
for n=1:N-1
    h = tspan(n+1)-tspan(n);
    k1 = odeFunc(t(n), y(:, n), Bmax, uL, uI, e, T, day, Ap);
    k2 = odeFunc(t(n) + 0.5*h, y(:,n) + 0.5.*k1*h, Bmax, uL, uI, e, T, day, Ap);
    k3 = odeFunc(t(n) + 0.5*h, y(:,n) + 0.5*k2*h, Bmax, uL, uI, e, T, day, Ap);
    k4 = odeFunc(t(n) + h, y(:,n) + k3*h, Bmax, uL, uI, e, T, day, Ap);
    y(:, n+1) = y(:,n)+((k1+(2*k2)+(2*k3)+k4))/6;
end
end

[CONNOR]

I apoligise I hadnt linked the files. Now the info file and the correct graph file have been linked

Answer 1

[David Goodmanson]

Hi Connor,

The reason your code is going to 10^26 is that you are missing an important factor on the last line of rk4, which is the width of the interval. You should have

y(:, n+1) = y(:,n)+((k1+(2*k2)+(2*k3)+k4))*h/6;

(running a test case on rk4 before trying to use it in the larger code might well have saved you time overall). After that change, you get a finite solution which takes care of the 10^26 effect.

The solution does not agree with the plot, and to agree with that plot you would have to adjust some of the starting values. Also it appears that the code builds up y0 as P,S,L,I,B but plots out y as S,L,I,R,P