I don't know why this code is saying, "Arrays have incompatible sizes for this operation."

This code has been giving me a lot of trouble trying to get it to work and now it is saying that the arrays are incompatable but I dont know why that is and where to look to fix that.

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; Ri; Bi];
[t,y] = rk4(@PSLIR, tspan, y0, Bmax, uL, uI, e, T, tspan, Ap);
%%      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(6);
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 = h*odeFunc(n, y(:, n), Bmax, uL, uI, e, T, day, Ap);
    k2 = h*odeFunc(n + 0.5, y(:,n) + 0.5.*k1, Bmax, uL, uI, e, T, day, Ap);
    k3 = h*odeFunc(n + 0.5, y(:,n) + 0.5*k2, Bmax, uL, uI, e, T, day, Ap);
    k4 = h*odeFunc(n + h, y(:,n) + k3, Bmax, uL, uI, e, T, day, Ap);
    y(:, n+1) = y(:,n)+((k1+(2*k2)+(2*k3)+k4))/6;
end
end

Answer 1

[Torsten]

Your vector of initial conditions is 6x1

y0 = [Pi; Si; Li; Ii; Ri; Bi];

but you only return a 5x1 vector in PSLIR

dydt = [dPdt; dSdt; dLdt; dIdt; dRdt];

Further, your first arguments in the calls to PSLIR are wrong (at least according to Runge-Kutta 4th order).

They must be t(n), t(n)+0.5*h, t(n)+0.5*h, t(n)+h in the calculation of k1,k2,k3 and k4 instead of n, n+0.5, n+0.5, n+h.

[CONNOR]

Thank you