MATLAB: Error: Unable to perform assignment because the left and right sides have a different number of elements.

arrayfor loop

The code shown below is my attempt to develop an SEIR model to solve the system of Ordinary Differential Equations using RK4. An error pops up when running the code stating :
"Unable to perform assignment because the left and right sides have a
different number of elements.
Error in SEIR_RK4 (line 55)
S(i+1) = S(i) + (s1 + 2*s2 + 2*s3 + s4)/6;"
I've looked through numerous sources but I'm still puzzeled how to make the code run. 
r0 = 3;
R = r0;%For Now
Tinc = 0.5; %Time of Incubation ~ 2 days
Tinf = 0.1; %Time of Infected ~ 10 days
t0 = 0; %Initial Time - Given
tf = 300;%Final Time - Given
N = 300; %Number of Points (or "Day Intervals") % M = N 
tol = 10^(-6); %Tolerance
h = (tf-t0)/N;%Step Size
t = [t0:h:tf]';%Discretized time step
S0 = 0.8; %Initial amount of Susceptibles
E0 = 0.2; %Initial amount of Incubants
I0 = 0.0; %Initial amount of Infectants
R0 = 0.0; %Initial amount of Recovered
%Initial Values X(0) - Matrix (4x1)
S = zeros(1,N); S(1) = S0;
E = zeros(1,N); E(1) = E0;
I = zeros(1,N); I(1) = I0;
R = zeros(1,N); R(1) = R0;
% f(t,x) = dx/dt
%f = @(t,x) [(-R*x(3)*x(1)/Tinf);(R*x(3)*x(1)/Tinf)-(x(2)/Tinc);(x(2)/Tinc)-(x(3)/Tinf);(x(3)/Tinf)];
dS = @(s,e,i,r) (-R.*i.*s/Tinf);
dE = @(s,e,i,r) (R.*i.*s/Tinf)-(e/Tinc);
dI = @(s,e,i,r) (e/Tinc)-(i/Tinf);
dR = @(s,e,i,r) (i/Tinf);
for i = 1:N
    %k1 = h*f(tn,yn)
    s1 = h*dS( S(i), E(i), I(i), R(i));
    e1 = h*dE( S(i), E(i), I(i), R(i));
    i1 = h*dI( S(i), E(i), I(i), R(i));
    r1 = h*dR( S(i), E(i), I(i), R(i));
    %k2 = h*f(tn + (h/2),yn + k1.*(h/2))
    s2 = h*dS( S(i) + 0.5*s1, E(i) + 0.5*e1, I(i) + 0.5*i1, R(i) + 0.5*r1);
    e2 = h*dE( S(i) + 0.5*s1, E(i) + 0.5*e1, I(i) + 0.5*i1, R(i) + 0.5*r1);
    i2 = h*dI( S(i) + 0.5*s1, E(i) + 0.5*e1, I(i) + 0.5*i1, R(i) + 0.5*r1);
    r2 = h*dR( S(i) + 0.5*s1, E(i) + 0.5*e1, I(i) + 0.5*i1, R(i) + 0.5*r1);
    %k3 = h*f(tn + (h/2),yn + k2.*(h/2))
    s3 = h*dS( S(i) + 0.5*s2, E(i) + 0.5*e2, I(i) + 0.5*i2, R(i) + 0.5*r2);
    e3 = h*dE( S(i) + 0.5*s2, E(i) + 0.5*e2, I(i) + 0.5*i2, R(i) + 0.5*r2);
    i3 = h*dI( S(i) + 0.5*s2, E(i) + 0.5*e2, I(i) + 0.5*i2, R(i) + 0.5*r2);
    r3 = h*dR( S(i) + 0.5*s2, E(i) + 0.5*e2, I(i) + 0.5*i2, R(i) + 0.5*r2);
    %k4 = h*f(tn + h,yn + h.*k3)
    s4 = h*dS( S(i) + s3, E(i) + e3, I(i) + i3, R(i) + r3);
    e4 = h*dE( S(i) + s3, E(i) + e3, I(i) + i3, R(i) + r3);
    i4 = h*dI( S(i) + s3, E(i) + e3, I(i) + i3, R(i) + r3);
    r4 = h*dR( S(i) + s3, E(i) + e3, I(i) + i3, R(i) + r3);
    %yn+1 = yn + (1/6)*h*(k1+2*k2+2*k3+k4)
    S(i+1) = S(i) + (s1 + 2*s2 + 2*s3 + s4)/6;
    E(i+1) = E(i) + (e1 + 2*e2 + 2*e3 + e4)/6;
    I(i+1) = I(i) + (i1 + 2*i2 + 2*i3 + i4)/6;
    R(i+1) = R(i) + (r1 + 2*r2 + 2*r3 + r4)/6;
end

Best Answer

Okay, so it turns out that you do not need to use cells. What you do need is to fix 
dS = @(s,e,i,r) (-R.*i.*s/Tinf);
R is not 
R = r0;%For Now
like you claim, because you later have
R = zeros(1,N); R(1) = R0;
so R is the vector for keeping track of the number of recovered individuals.
You need to rename one of the uses of R, as either the transmittability or the number of recovered individuals.
If you define
odefun = @(t,x) [dS(x(1),x(2),x(3),x(4)); dE(x(1),x(2),x(3),x(4)); dI(x(1),x(2),x(3),x(4)); dR(x(1),x(2),x(3),x(4))];
x0 = [S0; E0; I0; R0];
tspan = t;
[t, x] = ode23(odefun, tspan, x0);
then the results you get would suggest that your implementation of ode23 is incorrect. But also as you look at the E column (second column) and I column (third column) since the values cycle between negative and positive, that suggests that your dE and dI are incorrect.
Question: are your values intended to be proportions or intended to be absolute numbers? You do not have a population size, so it looks more like proportions ?
Related Question