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

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

[Daniel Pollard]

This code is extremely long, and as it's a function, I don't think I can easily run it as I don't know how to use it. Can you provide a minimum example which produces the same error?

I tried to make sense of it manually, but each variable is built from multiple others. "Different number of elements" means that lambda_S(j-1) has a different number of components to lambda_S(j)-(h/6)*(k1+2*k2+2*k3+k4);. I suspect the error is in the k's but I can't find it.

[Emman Obuya]

Dear Daniel unfortunately the code is a routine for backward-forward rK4 method which i got and fitted in my equations i may not be able to snapshot a smaller example but upon several checks i noticed it is only this part(lambda_S(j-1) = lambda_S(j)-(h/6)*(k1+2*k2+2*k3+k4);) that is bringing the error and i have failed to figure the way out..the reason am seeking for help

Answer 1

[Alan Stevens]

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

[Emman Obuya]

Dear Alan i looked into the comment you raised and made the definition, i was happy the comment i have been getting is nolonger appearing but then my code is not simply running could something else be a problem?? thanks for your reply in the meantime; below is the current code

function y =Try2( )
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=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);
    l10=min(l1,l2);
    l11=min(l3,l4);
    l12=min(l5,l6);
    l13=min(l7,l8);
    l14=min(l9,l10);  % l10  defined now
    l15=min(l11,l12);
    l16=min(l13,l14);
    l17=min(l15,l16); % defined now
    l18=min(temp1,l17);% "
    test=min(l15,l18); % "
    
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 ');

[Alan Stevens]

Your 'while' loop never ends! The variable 'test' just gets more and more negative. You need to decide if this is the right condition, or if 'test' is calculated correctly.

I put an iteration counter in this loop and got it to exit after a certain number of iterations. After you exit the loop, you come to the line

y(:,18)=E1;

Unfortunately, you haven't defined E1 anywhere.

Commenting this line out then allows the program to complete.

[Emman Obuya]

Dear Alan thanks again for the reply am really happy that am soon getting a way out from something that put me on hold for about 8 weeks now; sorry the "y(:,18)=E1;" was meant to be y(:,18)=V1; please am somehow new in matlab please do share with me the code that exited it will be of big help to me because i have no idea of how to put an iteration counter!!. Am just praying that the code exits! I have attached the correction of "E1"

function y =Try2( )
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=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);
    l10=min(l1,l2);
    l11=min(l3,l4);
    l12=min(l5,l6);
    l13=min(l7,l8);
    l14=min(l9,l10);  % l10 not defined
    l15=min(l11,l12);
    l16=min(l13,l14);
    l17=min(l15,l16); %not defined 
    l18=min(temp1,l17);% "
    test=min(l15,l18); % "
    
end
y(:,1)=t;
y(:,2)=S;
y(:,3)=V;
y(:,4)=I;
y(:,5)=R;
y(:,6)=B;
y(:,17)=S1;
y(:,18)=V1;
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 ');

[Alan Stevens]

Starting Immediately above the ‘while’ command, put

its = 0;

Then, change the ‘while’ command to, say,

while (test<0) && (its<100)

its = its + 1;

% ... rest of the code in the loop

end

However, you should still check your equations to see why test always seems to be negative.

Answer 2

[Emman Obuya]

Dear Alan, thanks for helping me out, checked the code and it only runs with the counter iteration you provided; Just have one last request; (U_note, U_one, U_two and U_three) are my 4 interventions from which am seeking an Optimal control but am seemingly not finding a code for plotting 2 situations(without and with control) in one plot say one of (S and S1; with say U_note, U_one while U_two=U_three=0 ) may you help me out with the part of the code? I believe one scenario would guide me thru other scenarios! thanks in advance

[Alan Stevens]

To plot say U_note and U_one on the same graph, do something like

plot(t,U_note,t,U_one), grid
xlabel('time'), ylabel('Unote and Uone')
legend('U_note','U_one')

Answer 3

[Emman Obuya]

Dear Alan the command

plot(t,S,t,S1), grid

xlabel('time'), ylabel('Susceptible')

legend('S','S1')

Surely plots the two lines in the same graph quite well, but the only part am missing now are the constraits like plotting (S, and S1) with (U_note, U_one and U_two) but U_three=0; Can i put those sevral constraints while plotting or they are supposed to be formed before i start plotting? Been thinking of a command like this;

plot(t,S,t,S1,'U_note','U_one','U_three'), grid

xlabel('time'), ylabel('Susceptible')

legend('S','S1')

But am getting an error

[Alan Stevens]

I'm not clear what you mean! You can plot other lines on the same graph by extending as follows:

plot(t,S,t,S1,t,U_note,t,U_one,t,U_three)

and then add labels in legend.

Try

doc plot

for more detail about using the plot command.

Answer 4

[Emman Obuya]

Am checking through; But as i look what i meant was that I am looking for plots with my interventions being (U_note=Vaccination, U_one= Education awareness, U_two= Treatment and U_three= Sanitation); I am trying to find plots (each with 2 lines of with and without control )of optimal control of (S,V,I R and B) for;

1. Figure 1; Plots of (SVIRB) when the interventions are binary for example(Vaccination+Treatment, or Sanitation+ Treatment).
2. Figure 2; Plots of (SVIRB) considering a combination of three intervention like (Vaccination+ Sanitation+ Treatment)
3. Figure 3, Plots of SVIRB with all the four interventions

It is this that am wondering whether i can put these scenarios in the command of plots since my code has all the situations, just wondering how to plot it out!!

Answer 5

[Emman Obuya]

Dear Alan, as a last question to my follow up, looking at the code i posted is there a way i can make these plots??

I am looking for plots with my interventions being (U_note=Vaccination, U_one= Education awareness, U_two= Treatment and U_three= Sanitation); I am trying to find plots (each with 2 lines of; with and without control )of optimal control of (S,V,I R and B) for;

1. Figure 1; Plots of (SVIRB) when the interventions are binary for example(Vaccination+Treatment, or Sanitation+ Treatment).
2. Figure 2; Plots of (SVIRB) considering a combination of three intervention like (Vaccination+ Sanitation+ Treatment)
3. Figure 3, Plots of SVIRB with all the four interventions

I beleive you might have plotted such figures?? please help me if you are able

[Alan Stevens]

I'm still not clear what you are trying to plot! The commands

plot(t,S,t,V,t,I,t,R,t,B),grid
xlabel('t'),ylabel('SVIRB')
legend('S','V','I','R','B')

produce the following:

That's after 100 iterations, but I can't tell if it's what you are looking for.

Answer 6

[Emman Obuya]

Thanks again Alan; Am almost there now but what i was meaning is a plot of this nature;

So am looking for plots of S, I, V, R and B with the two lines Red(without control) and Blue(With control) but when am varying the controls like for the first one like (Considering only 2 controls, then the next set considering like 3 controls being figure 2) !! in that order. I pray you really understand now

[Alan Stevens]

Well, these plot commands at the end of your program will produce those style lines:

figure
plot(t,S,'b',t,S1,'r'),grid
xlabel('Time');
ylabel('Susceptible');
legend('S','S1')
figure
plot(t,I,'b',t,I1,'r'),grid
xlabel('Time')
ylabel('Infected');
legend('I','I1')

However, they are a different shape from the ones you show above.