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

Question

Am finding the optimal control for the set of equations but am getting the same error in line 203

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

lambda_S(j-1) = lambda_S(j)-(h/6)*(k1+2*k2+2*k3+k4);

Have tried to check for similar answers but i seem not to get a way out easily. Your help will be of good benefit and thanks in advance

function y =Try( )
clc;
A1=0.02;
A2=10;
L1=1.6;
L2=4.0;
L3=0.6;
L4=0.8;
Capital_lambda=1000/100;
Alpha=0.5/100;
lambda=0.005/100;
beta=50/100;
k=100000;
b=0.0027/100;
d=1.5/100;
mu_h=0.006/100;
mu_b=2/100;
Gamma=10/100;
r=20/100;
test=-1;
N=10000;
t=linspace(0,100,N+1); % create a matrix starting from 0 going to 100 with 1001 values
h=1/N;
h2=h/2;
U_note=zeros(1,N+1); % Creating Zero matrices
U_one=zeros(1,N+1);
U_two=zeros(1,N+1);
U_three=zeros(1,N+1);
S=zeros(1,N+1);%S
S0=40;
S(1)=S0;
S1=zeros(1,N+1);%S1
S10=40;
S1(1)=S10;
V=zeros(1,N+1);%V
V0=30;
V(1)=V0;
V1=zeros(1,N+1);%V1
V10=30;
V1(1)=V10;
I=zeros(1,N+1);%I
I0=20;
I(1)=I0;
I1=zeros(1,N+1);%I1
I10=20;
I1(1)=I10;
R=zeros(1,N+1);%R

R0=10;
R(1)=R0;
R1=zeros(1,N+1);%R
R10=10;
R1(1)=R10;
B=zeros(1,N+1);%B

B0=20;
B(1)=B0;
B1=zeros(1,N+1);%B
B10=20;
B1(1)=B10;
lambda_S=zeros(1,N+1);
%lambda1(100)=0;
lambda_V=zeros(1,N+1);
%lambda2(100)=0;
lambda_I=zeros(1,N+1);
%lambda3(100)=0;
lambda_R=zeros(1,N+1);
%lambda4(100)=0;
lambda_B=zeros(1,N+1);
%lambda5(100)=0;
while(test<0) %Creates a loop to be repeated when the conditions are true
oldU_note=U_note;
oldU_one=U_one;
oldU_two=U_two;
oldU_three=U_three;
oldS=S;
oldS1=S1;
oldV=V;
oldV1=V1;
oldI=I;
oldI1=I1;
oldR=R;
oldR1=R1;
oldB=B;
oldB1=B1;
oldlambda_S=lambda_S;
oldlambda_V=lambda_V;
oldlambda_I=lambda_I;
oldlambda_R=lambda_R;
oldlambda_B=lambda_B;
for i=1:N
    k1=Capital_lambda+Gamma*R(i)+Alpha*V(i)-((1-U_one(i))*Alpha+mu_h+U_note(i))*S(i);
    k2=Capital_lambda+Gamma*(R(i)+h2*k1)+Alpha*(V(i)+h2*k1)-((1-0.5*(U_one(i)+U_one(i+1)))*Alpha+mu_h+0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k1));
    k3=Capital_lambda+Gamma*(R(i)+h2*k2)+Alpha*(V(i)+h2*k2)-((1-0.5*(U_one(i)+U_one(i+1)))*Alpha+mu_h+0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k2));
    k4=Capital_lambda+Gamma*(R(i)+h*k3)+Alpha*(V(i)+h*k3)-((1-U_one(i+1))*Alpha+mu_h+U_note(i+1))*(S(i)+h*k3);
    S(i+1)=S(i)+(h/6)*(k1+2*k2+2*k3+k4);
                                                                                                                  
 end
 for i=1:N
      k1=Capital_lambda+Gamma*R(i)+Alpha*V(i)-Alpha+mu_h;
      k2=Capital_lambda+Gamma*(R(i)+h2*k1)+Alpha*(V(i)+h2*k1)-Alpha+mu_h;
      k3=Capital_lambda+Gamma*(R(i)+h2*k2)+Alpha*(V(i)+h2*k2)-Alpha+mu_h;
      k4=Capital_lambda+Gamma*(R(i)+h*k3)+Alpha*(V(i)+h*k3)-Alpha+mu_h;
      S1(i+1)=S1(i)+(h/6)*(k1+2*k2+2*k3+k4);
                                                                                                                  
 end
 for i=1:N
     k1=U_note(i)*S(i)- mu_h*V(i)-Alpha*V(i);
     k2=0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k1)- mu_h*(V(i)+h2*k1)-Alpha*(V(i)+h2*k1);
     k3=0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k2)- mu_h*(V(i)+h2*k2)-Alpha*(V(i)+h2*k2);
     k4=0.5*(U_note(i+1))*(S(i)+h*k3)- mu_h*(V(i)+h*k3)-Alpha*(V(i)+h*k3);
     V(i+1)=V(i)+(h/6)*(k1+2*k2+2*k3+k4);
                                                                                                                 
 end   
 for i=1:N
     k1=-mu_h*V(i)-Alpha*V(i);
     k2=-mu_h*(V(i)+h2*k1)-Alpha*(V(i)+h2*k1);
     k3=-mu_h*(V(i)+h2*k2)-Alpha*(V(i)+h2*k2);
     k4=-mu_h*(V(i)+h*k3)-Alpha*(V(i)+h*k3);
     V1(i+1)=V1(i)+(h/6)*(k1+2*k2+2*k3+k4);
                                                                                                                  
 end
 for i=1:N
     k1=(1-U_one(i))*lambda*S(i)-((1-U_one(i))*beta + mu_h + d+ U_two(i)+r)*I(i);
     k2=(1-0.5*(U_one(i)+U_one(i+1)))*lambda*(S(i)+h2*k1)-((1-0.5*(U_one(i)+U_one(i+1)))*beta + mu_h + d+ 0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k1);
     k3=(1-0.5*(U_one(i)+U_one(i+1)))*lambda*(S(i)+h2*k2)-((1-0.5*(U_one(i)+U_one(i+1)))*beta + mu_h + d+ 0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k2);
     k4=(1-U_one(i+1))*lambda*(S(i)+h*k3)-((1-U_one(i+1))*beta + mu_h + d+U_two(i+1))+r*(I(i)+h*k3);
     I(i+1)=I(i)+(h/6)*(k1+2*k2+2*k3+k4);
                                                         
 end
                                                         
 for i=1:N
     k1=lambda*S(i)-beta + mu_h + d + r*I(i);
     k2=lambda*(S(i)+h2*k1)-beta + mu_h + d + r*(I(i)+h2*k1);
     k3=lambda*(S(i)+h2*k2)-beta + mu_h + d + r*(I(i)+h2*k2);
     k4=lambda*(S(i)+h*k3)-beta + mu_h + d + r*(I(i)+h*k3);
     I1(i+1)=I1(i)+(h/6)*(k1+2*k2+2*k3+k4);
 end
%fprintf('sum1=%f',sum1);
    for i=1:N
        k1=(U_two(i)+r)*I(i)-(mu_h+Gamma)*R(i);
        k2=(0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k1)-(mu_h+Gamma)*(R(i)+h2*k1);
        k3=(0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k2)-(mu_h+Gamma)*(R(i)+h2*k2);
        k4=(U_two(i+1))+r*(I(i)+h*k3)-(mu_h+Gamma)*(R(i)+h*k3);
        R(i+1)=R(i)+(h/6)*(k1+2*k2+2*k3+k4);
    end
    for i=1:N
        k1= r*I(i)-(mu_h+Gamma)*R(i);
        k2= r*(I(i)+h2*k1)-(mu_h+Gamma)*(R(i)+h2*k1);
        k3= r*(I(i)+h2*k2)-(mu_h+Gamma)*(R(i)+h2*k2);
        k4= r*(I(i)+h*k3)-(mu_h+Gamma)*(R(i)+h*k3);
        R1(i+1)=R1(i)+(h/6)*(k1+2*k2+2*k3+k4);
    end
    for i=1:N
        k1=beta*(1-U_one(i))*I(i)+((1-B(i)/k)*b*B(i))-(U_three(i)+mu_b)*B(i);
        k2=beta*((1-0.5*(U_one(i)+U_one(i+1)))*(I(i)+h2*k1)+((1-(B(i)+h2*k1)/k)*b*(B(i)+h2*k1))-(0.5*(U_three(i)+U_three(i+1))+mu_b)*(B(i)+h2*k1));
        k3=beta*((1-0.5*(U_one(i)+U_one(i+1)))*(I(i)+h2*k2)+((1-(B(i)+h2*k2)/k)*b*(B(i)+h2*k2))-(0.5*(U_three(i)+U_three(i+1))+mu_b)*(B(i)+h2*k2));
        k4=beta*((1-(U_one(i+1))*(I(i)+h*k3)+((1-(B(i)+h*k3)/k)*b*(B(i)+h*k3))-(U_three(i+1))+mu_b)*(B(i)+h*k3));
        B(i+1)=B(i)+(h/6)*(k1+2*k2+2*k3+k4);
    end
    for i=1:N
        k1=beta*I(i)+((1-B(i)/k)*b*B(i))-mu_b*B(i);
        k2=beta*(I(i)+h2*k1)+((1-(B(i)+h2*k1)/k)*b*(B(i)+h2*k1))-mu_b*(B(i)+h2*k1);
        k3=beta*(I(i)+h2*k2)+((1-(B(i)+h2*k2)/k)*b*(B(i)+h2*k2))-mu_b*(B(i)+h2*k2);
        k4=beta*(I(i)+h*k3)+((1-(B(i)+h*k3)/k)*b*(B(i)+h*k3))-mu_b*(B(i)+h*k3);
        B1(i+1)=B1(i)+(h/6)*(k1+2*k2+2*k3+k4);
    end
    for i=1:N
        j=N+2-i;
        k1=lambda_S(j)*((1-U_one(j))*lambda+mu_h+U_note(j))-A1-lambda_V*U_note(j)-lambda_I*lambda(1-U_one(j));
        k2=(lambda_S(j)-h2*k1)*((1-0.5*(U_one(j)+U_one(j-1)))*lambda+mu_h+0.5*(U_note(j)+U_note(j-1)))-A1-(lambda_V-h2*k1)*0.5*(U_note(j)+U_note(j-1))-(lambda_I-h2*k1)*lambda(1-0.5*(U_one(j)+U_one(j-1)));
        k3=(lambda_S(j)-h2*k2)*((1-0.5*(U_one(j)+U_one(j-1)))*lambda+mu_h+0.5*(U_note(j)+U_note(j-1)))-A1-(lambda_V-h2*k2)*0.5*(U_note(j)+U_note(j-1))-(lambda_I-h2*k2)*lambda(1-0.5*(U_one(j)+U_one(j-1)));
        k4=(lambda_S(j)-h*k3)*((1-U_one(j-1)))*lambda+mu_h+U_note(j-1)-A1-(lambda_V-h*k3)*U_note(j-1)-(lambda_I-h*k3)*lambda(1-U_one(j-1));
        lambda_S(j-1) = lambda_S(j)-(h/6)*(k1+2*k2+2*k3+k4);
        
    end
    
    for i=1:N
        j=N+2-i;
        k1=lambda_V(j)*(mu_h+Alpha)-lambda_S(j)*Alpha;
        k2=(lambda_V(j)-h2*k1)*(mu_h+Alpha)-(lambda_S(j)-h2*k1)*Alpha;
        k3=(lambda_V(j)-h2*k2)*(mu_h+Alpha)-(lambda_S(j)-h2*k2)*Alpha;
        k4=(lambda_V(j)-h*k3)*(mu_h+Alpha)-(lambda_S(j)-h*k3)*Alpha;
        lambda_V(j-1)=lambda_V(j)-(h/6)*(k1+2*k2+2*k3+k4);
    end
    
    
    for i=1:N
        j=N+2-i;
        k1=lambda_I(j)*(beta*(1-U_one(j))+mu_h+d+U_two(j))-A2-lambda_R(j)*Alpha*U_two(j)-lambda_B(j)*beta*(1-U_one(j));
        k2=(lambda_I(j)-h2*k1)*(beta*(1-0.5*(U_one(j)+U_one(j-1)))+mu_h+d+(U_two(j)+U_two(j-1)))-A2-(lambda_R(j)-h2*k1)*Alpha*(U_two(j)+U_two(j-1))-(lambda_B(j)-h2*k1)*beta*(1-0.5*(U_one(j)+U_one(j-1)));
        k3=(lambda_I(j)-h2*k2)*(beta*(1-0.5*(U_one(j)+U_one(j-1)))+mu_h+d+(U_two(j)+U_two(j-1)))-A2-(lambda_R(j)-h2*k2)*Alpha*(U_two(j)+U_two(j-1))-(lambda_B(j)-h2*k2)*beta*(1-0.5*(U_one(j)+U_one(j-1)));
        k4=(lambda_I(j)-h*k3)*(beta*(1-0.5*(U_one(j)+U_one(j-1)))+mu_h+d+(U_two(j)+U_two(j-1)))-A2-(lambda_R(j)-h*k3)*Alpha*(U_two(j)+U_two(j-1))-(lambda_B(j)-h*k3)*beta*(1-0.5*(U_one(j)+U_one(j-1)));
        lambda_I(j-1)=lambda_I(j)-(h/6)*(k1+2*k2(j)+2*k3(j)+k4(j));
    end
    
    for i=1:N
        j=N+2-i;
        k1=lambda_R(j)*(mu_h+Gamma)-lambda_S(j)*Gamma;
        k2=(lambda_R(j)-h2*k1)*(mu_h+Gamma)-(lambda_S(j)-h2*k1)*Gamma;
        k3=(lambda_R(j)-h2*k2)*(mu_h+Gamma)-(lambda_S(j)-h2*k2)*Gamma;
        k4=(lambda_R(j)-h*k3)*(mu_h+Gamma)-(lambda_S(j)-h*k3)*Gamma;
        lambda_R(j-1)=lambda_R(j)-(h/6)*(k1+2*k2+2*k3+k4);
    end
    for i=1:N
        j=N+2-i;
        k1=lambda_B*(mu_b+U_three);
        k2=(lambda_B(j)-h2*k1)*(mu_b+0.5*(U_three(j)+U_three(j-1)));
        k3=(lambda_B(j)-h2*k2)*(mu_b+0.5*(U_three(j)+U_three(j-1)));
        k4=(lambda_B(j)-h*k3)*(mu_b+0.5*(U_three(j)+U_three(j-1)));
        lambda_B(j-1)=lambda_B(j)-(h/6)*(k1+2*k2+2*k3+k4);
    end
    u0=((lambda_S-lambda_V)*S(i))/L1;
    U_note=0.5*(u0+oldU_note);
    u1=((lambda_I-lambda_S)*lambda*S(i)-lambda_I*beta)/L2;
    U_one=0.5*(u1+oldU_one);
    u2=((lambda_I-lambda_R)*I(i))/L3;
    U_two=0.5*(u2+oldU_two);
    u3=(lambda_B*B(i))/L4;
    U_three=0.5*(u3+oldU_three);
    temp1=0.001*sum(abs(U_note))-sum(abs(oldU_note-U_note));
    temp2=0.001*sum(abs(U_one))-sum(abs(oldU_one-U_one));
    temp3=0.001*sum(abs(U_two))-sum(abs(oldU_two-U_two));
    temp4=0.001*sum(abs(U_three))-sum(abs(oldU_three-U_three));
    temp5=0.001*sum(abs(S))-sum(abs(oldS-S));
    temp6=0.001*sum(abs(V))-sum(abs(oldV-V));
    temp7=0.001*sum(abs(I))-sum(abs(oldI-I));
    temp8=0.001*sum(abs(R))-sum(abs(oldR-R));
    temp9=0.001*sum(abs(B))-sum(abs(oldB-B));
    temp10=0.001*sum(abs(lambda_S))-sum(abs(oldlambda_S-lambda_S));
    temp11=0.001*sum(abs(lambda_V))-sum(abs(oldlambda_V-lambda_V));
    temp12=0.001*sum(abs(lambda_I))-sum(abs(oldlambda_I-lambda_I));
    temp13=0.001*sum(abs(lambda_R))-sum(abs(oldlambda_R-lambda_R));
    temp14=0.001*sum(abs(lambda_B))-sum(abs(oldlambda_B-lambda_B));
    temp15=0.001*sum(abs(S1))-sum(abs(oldS1-S1));
    temp16=0.001*sum(abs(V1))-sum(abs(oldV1-V1));
    temp17=0.001*sum(abs(I1))-sum(abs(oldI1-I1));
    temp18=0.001*sum(abs(R1))-sum(abs(oldR1-R1));
    temp19=0.001*sum(abs(B1))-sum(abs(oldB1-B1));
    l1=min(temp2,temp3);
    l2=min(temp4,temp5);
    l3=min(temp6,temp7);
    l4=min(temp8,temp9);
    l5=min(temp10,temp11);
    l6=min(temp12,temp13);
    l7=min(temp14,temp15);
    l8=min(temp16,temp17);
    l9=min(temp18,temp19);
    l11=min(l1,l2);
    l12=min(l3,l4);
    l13=min(l5,l6);
    l14=min(l7,l8);
    l15=min(l9,l10);
    l16=min(l11,l12);
    l17=min(l13,l14);
    l18=min(l15,l16);
    l19=min(l17,l18);
    l20=min(temp1,l19);
    test=min(l15,l20);
    
end
y(:,1)=t;
y(:,2)=S;
y(:,3)=V;
y(:,4)=I;
y(:,5)=R;
y(:,6)=B;
y(:,17)=S1;
y(:,18)=E1;
y(:,19)=I1;
y(:,20)=R1;
y(:,21)=B1;
y(:,14)=U_note;
y(:,15)=U_one;
y(:,16)=U_two;
y(:,17)=U_three;
y(:,8)=lambda_S;
y(:,9)=lambda_V;
y(:,10)=lambda_I;
y(:,11)=lambda_R;
y(:,12)=lambda_B;
plot(t,S);
xlabel('Time');
ylabel('Susceptible human');
plot(t,V);
xlabel('Time')
ylabel('Exposed human');
plot(t,R);
xlabel('Time')
ylabel('Exposed ');

Answer 1

There are a lot of missing indices and definitions. In the following the indices are repaired; however, you will need to define l10 (see line 256 or thereabouts) to get the routine to work

A1=0.02;
A2=10;
L1=1.6;
L2=4.0;
L3=0.6;
L4=0.8;
Capital_lambda=1000/100;
Alpha=0.5/100;
lambda=0.005/100;
beta=50/100;
k=100000;
b=0.0027/100;
d=1.5/100;
mu_h=0.006/100;
mu_b=2/100;
Gamma=10/100;
r=20/100;
test=-1;
N=1000;
t=linspace(0,100,N+1); % create a matrix starting from 0 going to 100 with 1001 values
h=1/N;
h2=h/2;
U_note=zeros(1,N+1); % Creating Zero matrices
U_one=zeros(1,N+1);
U_two=zeros(1,N+1);
U_three=zeros(1,N+1);
S=zeros(1,N+1);%S
S0=40;
S(1)=S0;
S1=zeros(1,N+1);%S1
S10=40;
S1(1)=S10;
V=zeros(1,N+1);%V
V0=30;
V(1)=V0;
V1=zeros(1,N+1);%V1
V10=30;
V1(1)=V10;
I=zeros(1,N+1);%I
I0=20;
I(1)=I0;
I1=zeros(1,N+1);%I1
I10=20;
I1(1)=I10;
R=zeros(1,N+1);%R

R0=10;
R(1)=R0;
R1=zeros(1,N+1);%R
R10=10;
R1(1)=R10;
B=zeros(1,N+1);%B

B0=20;
B(1)=B0;
B1=zeros(1,N+1);%B
B10=20;
B1(1)=B10;
lambda_S=zeros(1,N+1);
%lambda1(100)=0;
lambda_V=zeros(1,N+1);
%lambda2(100)=0;
lambda_I=zeros(1,N+1);
%lambda3(100)=0;
lambda_R=zeros(1,N+1);
%lambda4(100)=0;
lambda_B=zeros(1,N+1);
%lambda5(100)=0;
while(test<0) %Creates a loop to be repeated when the conditions are true
oldU_note=U_note;
oldU_one=U_one;
oldU_two=U_two;
oldU_three=U_three;
oldS=S;
oldS1=S1;
oldV=V;
oldV1=V1;
oldI=I;
oldI1=I1;
oldR=R;
oldR1=R1;
oldB=B;
oldB1=B1;
oldlambda_S=lambda_S;
oldlambda_V=lambda_V;
oldlambda_I=lambda_I;
oldlambda_R=lambda_R;
oldlambda_B=lambda_B;
for i=1:N
    k1=Capital_lambda+Gamma*R(i)+Alpha*V(i)-((1-U_one(i))*Alpha+mu_h+U_note(i))*S(i);
    k2=Capital_lambda+Gamma*(R(i)+h2*k1)+Alpha*(V(i)+h2*k1)-((1-0.5*(U_one(i)+U_one(i+1)))*Alpha+mu_h+0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k1));
    k3=Capital_lambda+Gamma*(R(i)+h2*k2)+Alpha*(V(i)+h2*k2)-((1-0.5*(U_one(i)+U_one(i+1)))*Alpha+mu_h+0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k2));
    k4=Capital_lambda+Gamma*(R(i)+h*k3)+Alpha*(V(i)+h*k3)-((1-U_one(i+1))*Alpha+mu_h+U_note(i+1))*(S(i)+h*k3);
    S(i+1)=S(i)+(h/6)*(k1+2*k2+2*k3+k4);
                                                                                                                  
 end
 for i=1:N
      k1=Capital_lambda+Gamma*R(i)+Alpha*V(i)-Alpha+mu_h;
      k2=Capital_lambda+Gamma*(R(i)+h2*k1)+Alpha*(V(i)+h2*k1)-Alpha+mu_h;
      k3=Capital_lambda+Gamma*(R(i)+h2*k2)+Alpha*(V(i)+h2*k2)-Alpha+mu_h;
      k4=Capital_lambda+Gamma*(R(i)+h*k3)+Alpha*(V(i)+h*k3)-Alpha+mu_h;
      S1(i+1)=S1(i)+(h/6)*(k1+2*k2+2*k3+k4);
                                                                                                                  
 end
 for i=1:N
     k1=U_note(i)*S(i)- mu_h*V(i)-Alpha*V(i);
     k2=0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k1)- mu_h*(V(i)+h2*k1)-Alpha*(V(i)+h2*k1);
     k3=0.5*(U_note(i)+U_note(i+1))*(S(i)+h2*k2)- mu_h*(V(i)+h2*k2)-Alpha*(V(i)+h2*k2);
     k4=0.5*(U_note(i+1))*(S(i)+h*k3)- mu_h*(V(i)+h*k3)-Alpha*(V(i)+h*k3);
     V(i+1)=V(i)+(h/6)*(k1+2*k2+2*k3+k4);
                                                                                                                 
 end   
 for i=1:N
     k1=-mu_h*V(i)-Alpha*V(i);
     k2=-mu_h*(V(i)+h2*k1)-Alpha*(V(i)+h2*k1);
     k3=-mu_h*(V(i)+h2*k2)-Alpha*(V(i)+h2*k2);
     k4=-mu_h*(V(i)+h*k3)-Alpha*(V(i)+h*k3);
     V1(i+1)=V1(i)+(h/6)*(k1+2*k2+2*k3+k4);
                                                                                                                  
 end
 for i=1:N
     k1=(1-U_one(i))*lambda*S(i)-((1-U_one(i))*beta + mu_h + d+ U_two(i)+r)*I(i);
     k2=(1-0.5*(U_one(i)+U_one(i+1)))*lambda*(S(i)+h2*k1)-((1-0.5*(U_one(i)+U_one(i+1)))*beta + mu_h + d+ 0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k1);
     k3=(1-0.5*(U_one(i)+U_one(i+1)))*lambda*(S(i)+h2*k2)-((1-0.5*(U_one(i)+U_one(i+1)))*beta + mu_h + d+ 0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k2);
     k4=(1-U_one(i+1))*lambda*(S(i)+h*k3)-((1-U_one(i+1))*beta + mu_h + d+U_two(i+1))+r*(I(i)+h*k3);
     I(i+1)=I(i)+(h/6)*(k1+2*k2+2*k3+k4);
                                                         
 end
                                                         
 for i=1:N
     k1=lambda*S(i)-beta + mu_h + d + r*I(i);
     k2=lambda*(S(i)+h2*k1)-beta + mu_h + d + r*(I(i)+h2*k1);
     k3=lambda*(S(i)+h2*k2)-beta + mu_h + d + r*(I(i)+h2*k2);
     k4=lambda*(S(i)+h*k3)-beta + mu_h + d + r*(I(i)+h*k3);
     I1(i+1)=I1(i)+(h/6)*(k1+2*k2+2*k3+k4);
 end
%fprintf('sum1=%f',sum1);
    for i=1:N
        k1=(U_two(i)+r)*I(i)-(mu_h+Gamma)*R(i);
        k2=(0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k1)-(mu_h+Gamma)*(R(i)+h2*k1);
        k3=(0.5*(U_two(i)+U_two(i+1))+r)*(I(i)+h2*k2)-(mu_h+Gamma)*(R(i)+h2*k2);
        k4=(U_two(i+1))+r*(I(i)+h*k3)-(mu_h+Gamma)*(R(i)+h*k3);
        R(i+1)=R(i)+(h/6)*(k1+2*k2+2*k3+k4);
    end
    for i=1:N
        k1= r*I(i)-(mu_h+Gamma)*R(i);
        k2= r*(I(i)+h2*k1)-(mu_h+Gamma)*(R(i)+h2*k1);
        k3= r*(I(i)+h2*k2)-(mu_h+Gamma)*(R(i)+h2*k2);
        k4= r*(I(i)+h*k3)-(mu_h+Gamma)*(R(i)+h*k3);
        R1(i+1)=R1(i)+(h/6)*(k1+2*k2+2*k3+k4);
    end
    for i=1:N
        k1=beta*(1-U_one(i))*I(i)+((1-B(i)/k)*b*B(i))-(U_three(i)+mu_b)*B(i);
        k2=beta*((1-0.5*(U_one(i)+U_one(i+1)))*(I(i)+h2*k1)+((1-(B(i)+h2*k1)/k)*b*(B(i)+h2*k1))-(0.5*(U_three(i)+U_three(i+1))+mu_b)*(B(i)+h2*k1));
        k3=beta*((1-0.5*(U_one(i)+U_one(i+1)))*(I(i)+h2*k2)+((1-(B(i)+h2*k2)/k)*b*(B(i)+h2*k2))-(0.5*(U_three(i)+U_three(i+1))+mu_b)*(B(i)+h2*k2));
        k4=beta*((1-(U_one(i+1))*(I(i)+h*k3)+((1-(B(i)+h*k3)/k)*b*(B(i)+h*k3))-(U_three(i+1))+mu_b)*(B(i)+h*k3));
        B(i+1)=B(i)+(h/6)*(k1+2*k2+2*k3+k4);
    end
    for i=1:N
        k1=beta*I(i)+((1-B(i)/k)*b*B(i))-mu_b*B(i);
        k2=beta*(I(i)+h2*k1)+((1-(B(i)+h2*k1)/k)*b*(B(i)+h2*k1))-mu_b*(B(i)+h2*k1);
        k3=beta*(I(i)+h2*k2)+((1-(B(i)+h2*k2)/k)*b*(B(i)+h2*k2))-mu_b*(B(i)+h2*k2);
        k4=beta*(I(i)+h*k3)+((1-(B(i)+h*k3)/k)*b*(B(i)+h*k3))-mu_b*(B(i)+h*k3);
        B1(i+1)=B1(i)+(h/6)*(k1+2*k2+2*k3+k4);
    end
    for i=1:N
        j=N+2-i;
        k1=lambda_S(j)*((1-U_one(j))*lambda+mu_h+U_note(j))-A1-lambda_V(j)*U_note(j)-lambda_I(j)*lambda*(1-U_one(j));
        k2=(lambda_S(j)-h2*k1)*((1-0.5*(U_one(j)+U_one(j-1)))*lambda+mu_h+0.5*(U_note(j)+U_note(j-1)))-A1-(lambda_V(j)-h2*k1)*0.5*(U_note(j)+U_note(j-1))-(lambda_I(j)-h2*k1)*lambda*(1-0.5*(U_one(j)+U_one(j-1)));
        k3=(lambda_S(j)-h2*k2)*((1-0.5*(U_one(j)+U_one(j-1)))*lambda+mu_h+0.5*(U_note(j)+U_note(j-1)))-A1-(lambda_V(j)-h2*k2)*0.5*(U_note(j)+U_note(j-1))-(lambda_I(j)-h2*k2)*lambda*(1-0.5*(U_one(j)+U_one(j-1)));
        k4=(lambda_S(j)-h*k3)*((1-U_one(j-1)))*lambda+mu_h+U_note(j-1)-A1-(lambda_V(j)-h*k3)*U_note(j-1)-(lambda_I(j)-h*k3)*lambda*(1-U_one(j-1));
        lambda_S(j-1) = lambda_S(j)-(h/6)*(k1+2*k2+2*k3+k4);
        
    end
    
    for i=1:N
        j=N+2-i;
        k1=lambda_V(j)*(mu_h+Alpha)-lambda_S(j)*Alpha;
        k2=(lambda_V(j)-h2*k1)*(mu_h+Alpha)-(lambda_S(j)-h2*k1)*Alpha;
        k3=(lambda_V(j)-h2*k2)*(mu_h+Alpha)-(lambda_S(j)-h2*k2)*Alpha;
        k4=(lambda_V(j)-h*k3)*(mu_h+Alpha)-(lambda_S(j)-h*k3)*Alpha;
        lambda_V(j-1)=lambda_V(j)-(h/6)*(k1+2*k2+2*k3+k4);
    end
    
    
    for i=1:N
        j=N+2-i;
        k1=lambda_I(j)*(beta*(1-U_one(j))+mu_h+d+U_two(j))-A2-lambda_R(j)*Alpha*U_two(j)-lambda_B(j)*beta*(1-U_one(j));
        k2=(lambda_I(j)-h2*k1)*(beta*(1-0.5*(U_one(j)+U_one(j-1)))+mu_h+d+(U_two(j)+U_two(j-1)))-A2-(lambda_R(j)-h2*k1)*Alpha*(U_two(j)+U_two(j-1))-(lambda_B(j)-h2*k1)*beta*(1-0.5*(U_one(j)+U_one(j-1)));
        k3=(lambda_I(j)-h2*k2)*(beta*(1-0.5*(U_one(j)+U_one(j-1)))+mu_h+d+(U_two(j)+U_two(j-1)))-A2-(lambda_R(j)-h2*k2)*Alpha*(U_two(j)+U_two(j-1))-(lambda_B(j)-h2*k2)*beta*(1-0.5*(U_one(j)+U_one(j-1)));
        k4=(lambda_I(j)-h*k3)*(beta*(1-0.5*(U_one(j)+U_one(j-1)))+mu_h+d+(U_two(j)+U_two(j-1)))-A2-(lambda_R(j)-h*k3)*Alpha*(U_two(j)+U_two(j-1))-(lambda_B(j)-h*k3)*beta*(1-0.5*(U_one(j)+U_one(j-1)));
        lambda_I(j-1)=lambda_I(j)-(h/6)*(k1+2*k2+2*k3+k4);
    end
    
    for i=1:N
        j=N+2-i;
        k1=lambda_R(j)*(mu_h+Gamma)-lambda_S(j)*Gamma;
        k2=(lambda_R(j)-h2*k1)*(mu_h+Gamma)-(lambda_S(j)-h2*k1)*Gamma;
        k3=(lambda_R(j)-h2*k2)*(mu_h+Gamma)-(lambda_S(j)-h2*k2)*Gamma;
        k4=(lambda_R(j)-h*k3)*(mu_h+Gamma)-(lambda_S(j)-h*k3)*Gamma;
        lambda_R(j-1)=lambda_R(j)-(h/6)*(k1+2*k2+2*k3+k4);
    end
    for i=1:N
        j=N+2-i;
        k1=lambda_B(j)*(mu_b+U_three(j));
        k2=(lambda_B(j)-h2*k1)*(mu_b+0.5*(U_three(j)+U_three(j-1)));
        k3=(lambda_B(j)-h2*k2)*(mu_b+0.5*(U_three(j)+U_three(j-1)));
        k4=(lambda_B(j)-h*k3)*(mu_b+0.5*(U_three(j)+U_three(j-1)));
        lambda_B(j-1)=lambda_B(j)-(h/6)*(k1+2*k2+2*k3+k4);
    end
    u0=((lambda_S-lambda_V)*S(i))/L1;
    U_note=0.5*(u0+oldU_note);
    u1=((lambda_I-lambda_S)*lambda*S(i)-lambda_I*beta)/L2;
    U_one=0.5*(u1+oldU_one);
    u2=((lambda_I-lambda_R)*I(i))/L3;
    U_two=0.5*(u2+oldU_two);
    u3=(lambda_B*B(i))/L4;
    U_three=0.5*(u3+oldU_three);
    temp1=0.001*sum(abs(U_note))-sum(abs(oldU_note-U_note));
    temp2=0.001*sum(abs(U_one))-sum(abs(oldU_one-U_one));
    temp3=0.001*sum(abs(U_two))-sum(abs(oldU_two-U_two));
    temp4=0.001*sum(abs(U_three))-sum(abs(oldU_three-U_three));
    temp5=0.001*sum(abs(S))-sum(abs(oldS-S));
    temp6=0.001*sum(abs(V))-sum(abs(oldV-V));
    temp7=0.001*sum(abs(I))-sum(abs(oldI-I));
    temp8=0.001*sum(abs(R))-sum(abs(oldR-R));
    temp9=0.001*sum(abs(B))-sum(abs(oldB-B));
    temp10=0.001*sum(abs(lambda_S))-sum(abs(oldlambda_S-lambda_S));
    temp11=0.001*sum(abs(lambda_V))-sum(abs(oldlambda_V-lambda_V));
    temp12=0.001*sum(abs(lambda_I))-sum(abs(oldlambda_I-lambda_I));
    temp13=0.001*sum(abs(lambda_R))-sum(abs(oldlambda_R-lambda_R));
    temp14=0.001*sum(abs(lambda_B))-sum(abs(oldlambda_B-lambda_B));
    temp15=0.001*sum(abs(S1))-sum(abs(oldS1-S1));
    temp16=0.001*sum(abs(V1))-sum(abs(oldV1-V1));
    temp17=0.001*sum(abs(I1))-sum(abs(oldI1-I1));
    temp18=0.001*sum(abs(R1))-sum(abs(oldR1-R1));
    temp19=0.001*sum(abs(B1))-sum(abs(oldB1-B1));
    l1=min(temp2,temp3);
    l2=min(temp4,temp5);
    l3=min(temp6,temp7);
    l4=min(temp8,temp9);
    l5=min(temp10,temp11);
    l6=min(temp12,temp13);
    l7=min(temp14,temp15);
    l8=min(temp16,temp17);
    l9=min(temp18,temp19);
    l11=min(l1,l2);
    l12=min(l3,l4);
    l13=min(l5,l6);
    l14=min(l7,l8);
    l15=min(l9,l10);  % l10 not defined
    l16=min(l11,l12);
    l17=min(l13,l14);
    l18=min(l15,l16); %not defined 
    l19=min(l17,l18); %  "
    l20=min(temp1,l19);% "

    test=min(l15,l20); % "
    
end
y(:,1)=t;
y(:,2)=S;
y(:,3)=V;
y(:,4)=I;
y(:,5)=R;
y(:,6)=B;
y(:,17)=S1;
y(:,18)=E1;
y(:,19)=I1;
y(:,20)=R1;
y(:,21)=B1;
y(:,14)=U_note;
y(:,15)=U_one;
y(:,16)=U_two;
y(:,17)=U_three;
y(:,8)=lambda_S;
y(:,9)=lambda_V;
y(:,10)=lambda_I;
y(:,11)=lambda_R;
y(:,12)=lambda_B;
plot(t,S);
xlabel('Time');
ylabel('Susceptible human');
plot(t,V);
xlabel('Time')
ylabel('Exposed human');
plot(t,R);
xlabel('Time')
ylabel('Exposed ');