test = -1;
tf=10;
delta = 0.001;
M = 50;
t = linspace(0,tf,M+1);
h = tf/M;
h2 = h/2;
pi=0.99999999999999999999;
beta1=0.112;
beta2=0.092;
omega=0.121;
mo=0.43;
m1=0.23;
mu=0.421;
beta3=0.054;
Ko=0.311;
d1=0.15;
eta1=0.33;
alpha1=0.22;
d2=0.17;
eta2=0.30;
alpha2=0.19;
delta1=0.002;
G0=0.00123;
G1=0.00145;
G2=0.00211;
K1=0.23;
K2=0.01;
K3=0.21;
x1=zeros(1,M+1);
x2=zeros(1,M+1);
x3=zeros(1,M+1);
x4=zeros(1,M+1);
x5=zeros(1,M+1);
u1 = zeros(1,M+1); u2 = zeros(1,M+1); u3 = zeros(1,M+1);
lambda1 = zeros(1,M+1);
lambda2 = zeros(1,M+1);
lambda3 = zeros(1,M+1);
lambda4 = zeros(1,M+1);
lambda5 = zeros(1,M+1);
x1(1) = 20;
x2(1) = 20;
x3(1) = 20;
x4(1) = 20;
x5(1) = 5;
while(test < 0)
oldu1 = u1;
oldu2 = u2;
oldu3 = u3;
oldx1 = x1;
oldx2 = x2;
oldx3 = x3;
oldx4 = x4;
oldx5 = x5;
oldlambda1 = lambda1;
oldlambda2 = lambda2;
oldlambda3 = lambda3;
oldlambda4 = lambda4;
oldlambda5 = lambda5;
for i = 1:M
m11 = pi-((1-u1(i)).*beta1.*x1(i).*x2(i))/(1+mo.*x2(i))-((1-u2(i)).*beta2.*x1(i).*x3(i))/(1+m1.*x3(i))-mu.*x1(i)-((1-u3(i)).*beta3.*x1(i).*x4(i))/(x4(i)+Ko);
m12 = ((1-u1(i)).*beta1.*x1(i).*x2(i))/(1+mo.*x2(i))-(d1+u2(i).*eta1+u3(i).*alpha1).*x2(i);
m13 = ((1-u2(i)).*beta2.*x1(i).*x3(i))/(1+m1.*x3(i))-(d2+u2(i).*eta2+u3(i).*alpha2).*x3(i);
m14 = u3(i).*alpha1.*x4(i)+u3(i).*alpha2.*x3(i)-delta1.*x4(i);
m15 = u2(i).*eta1.*x2(i)+u2(i).*eta2.*x3(i)-mu.*x5(i);
m21 = pi-((1-(0.5.*(u1(i)+u1(i+1)).*beta1.*(x1(i)+h2.*m11).*(x2(i)+h2.*m12))))/(1+mo.*(x2(i)+h2.*m12))-((1-(u2(i)+u2(i+1)).*beta2.*(x1(i)+h2.*m11).*(x3(i)+h2.*m13)))/(1+m1.*(x3(i)+h2.*m13))-mu.*(x1(i)+h2.*m11)-((1-(u3(i)+u3(i+1)).*beta3.*(x1(i)+h2.*m11).*(x4(i)+h2.*m14)))/((x4(i)+h2.*m14)+Ko);
m22 = ((1-(0.5.*(u1(i)+u1(i+1)).*beta1.*(x1(i)+h2.*m11).*(x2(i)+h2.*m12))))/(1+mo.*(x2(i)+h2.*m12))-(d1+(u2(i)+u2(i+1)).*eta1+(u3(i)+u3(i+1)).*alpha1).*(x2(i)+h2.*m12);
m23 = ((1-(u2(i)+u2(i+1)).*beta2.*(x1(i)+h2.*m11).*(x3(i)+h2.*m13)))/(1+m1.*(x3(i)+h2.*m13))-(d2+(u2(i)+u2(i+1)).*eta2+(u3(i)+u3(i+1)).*alpha2).*(x3(i)+h2.*m13);
m24 = (u3(i)+u3(i+1)).*alpha1.*(x4(i)+h2.*m14)+(u3(i)+u3(i+1)).*alpha2.*(x3(i)+h2.*m13)-delta1.*(x4(i)+h2.*m14);
m25 = (u2(i)+u2(i+1)).*eta1.*(x2(i)+h2.*m12)+(u2(i)+u2(i+1)).*eta2.*(x3(i)+h2.*m13)-mu.*(x5(i)+h2.*m15);
m31 = pi-((1-(0.5.*(u1(i)+u1(i+1)).*beta1.*(x1(i)+h2.*m21).*(x2(i)+h2.*m22))))/(1+mo.*(x2(i)+h2.*m22))-((1-(u2(i)+u2(i+1)).*beta2.*(x1(i)+h2.*m21).*(x3(i)+h2.*m23)))/(1+m1.*(x3(i)+h2.*m23))-mu.*(x1(i)+h2.*m21)-((1-(u3(i)+u3(i+1)).*beta3.*(x1(i)+h2.*m21).*(x4(i)+h2.*m24)))/((x4(i)+h2.*m24)+Ko);
m32 = ((1-(0.5.*(u1(i)+u1(i+1)).*beta1.*(x1(i)+h2.*m21).*(x2(i)+h2.*m22))))/(1+mo.*(x2(i)+h2.*m22))-(d1+(u2(i)+u2(i+1)).*eta1+(u3(i)+u3(i+1)).*alpha1).*(x2(i)+h2.*m22);
m33 = ((1-(u2(i)+u2(i+1)).*beta2.*(x1(i)+h2.*m21).*(x3(i)+h2.*m23)))/(1+m1.*(x3(i)+h2.*m23))-(d2+(u2(i)+u2(i+1)).*eta2+(u3(i)+u3(i+1)).*alpha2).*(x3(i)+h2.*m23);
m34 = (u3(i)+u3(i+1)).*alpha1.*(x4(i)+h2.*m24)+(u3(i)+u3(i+1)).*alpha2.*(x3(i)+h2.*m23)-delta1.*(x4(i)+h2.*m24);
m35 = (u2(i)+u2(i+1)).*eta1.*(x2(i)+h2.*m22)+(u2(i)+u2(i+1)).*eta2.*(x3(i)+h2.*m23)-mu.*(x5(i)+h2.*m25);
m41 = pi-((1-(0.5.*(u1(i)+u1(i+1)).*beta1.*(x1(i)+h2.*m31).*(x2(i)+h2.*m32))))/(1+mo.*(x2(i)+h2.*m32))-((1-(u2(i)+u2(i+1)).*beta2.*(x1(i)+h2.*m31).*(x3(i)+h2.*m33)))/(1+m1.*(x3(i)+h2.*m33))-mu.*(x1(i)+h2.*m31)-((1-(u3(i)+u3(i+1)).*beta3.*(x1(i)+h2.*m31).*(x4(i)+h2.*m34)))/((x4(i)+h2.*m34)+Ko);
m42 = ((1-(0.5.*(u1(i)+u1(i+1)).*beta1.*(x1(i)+h2.*m31).*(x2(i)+h2.*m32))))/(1+mo.*(x2(i)+h2.*m32))-(d1+(u2(i)+u2(i+1)).*eta1+(u3(i)+u3(i+1)).*alpha1).*(x2(i)+h2.*m32);
m43 = ((1-(u2(i)+u2(i+1)).*beta2.*(x1(i)+h2.*m31).*(x3(i)+h2.*m33)))/(1+m1.*(x3(i)+h2.*m33))-(d2+(u2(i)+u2(i+1)).*eta2+(u3(i)+u3(i+1)).*alpha2).*(x3(i)+h2.*m33);
m44 = (u3(i)+u3(i+1)).*alpha1.*(x4(i)+h2.*m34)+(u3(i)+u3(i+1)).*alpha2.*(x3(i)+h2.*m33)-delta1.*(x4(i)+h2.*m34);
m45 = (u2(i)+u2(i+1)).*eta1.*(x2(i)+h2.*m32)+(u2(i)+u2(i+1)).*eta2.*(x3(i)+h2.*m33)-mu.*(x5(i)+h2.*m35);
m51 = pi-((1-(0.5.*(u1(i)+u1(i+1)).*beta1.*(x1(i)+h2.*m41).*(x2(i)+h2.*m42))))/(1+mo.*(x2(i)+h2.*m42))-((1-(u2(i)+u2(i+1)).*beta2.*(x1(i)+h2.*m41).*(x3(i)+h2.*m43)))/(1+m1.*(x3(i)+h2.*m43))-mu.*(x1(i)+h2.*m41)-((1-(u3(i)+u3(i+1)).*beta3.*(x1(i)+h2.*m41).*(x4(i)+h2.*m44)))/((x4(i)+h2.*m44)+Ko);
m52 = ((1-(0.5.*(u1(i)+u1(i+1)).*beta1.*(x1(i)+h2.*m41).*(x2(i)+h2.*m42))))/(1+mo.*(x2(i)+h2.*m42))-(d1+(u2(i)+u2(i+1)).*eta1+(u3(i)+u3(i+1)).*alpha1).*(x2(i)+h2.*m42);
m53 = ((1-(u2(i)+u2(i+1)).*beta2.*(x1(i)+h2.*m41).*(x3(i)+h2.*m43)))/(1+m1.*(x3(i)+h2.*m43))-(d2+(u2(i)+u2(i+1)).*eta2+(u3(i)+u3(i+1)).*alpha2).*(x3(i)+h2.*m43);
m54 = (u3(i)+u3(i+1)).*alpha1.*(x4(i)+h2.*m44)+(u3(i)+u3(i+1)).*alpha2.*(x3(i)+h2.*m43)-delta1.*(x4(i)+h2.*m44);
m55 = (u2(i)+u2(i+1)).*eta1.*(x2(i)+h2.*m42)+(u2(i)+u2(i+1)).*eta2.*(x3(i)+h2.*m43)-mu.*(x5(i)+h2.*m45);
x1(i+1) = x1(i) + (h/6)*(m11 + 2*m21 + 2*m31 + 2.*m41 + m51);
x2(i+1) = x2(i) + (h/6)*(m12 + 2*m22 + 2*m32 + 2.*m42 + m52);
x3(i+1) = x3(i) + (h/6)*(m13 + 2*m23 + 2*m33 + 2.*m43 + m53);
x4(i+1) = x4(i) + (h/6)*(m14 + 2*m24 + 2*m34 + 2.*m44 + m54);
x5(i+1) = x5(i) + (h/6)*(m15 + 2.*m25 + 2.*m35 + 2.*m45 + m55);
end
for i=1:M;
j=M+2-1;
n11=(lambda2 - lambda1).*((1-u1(j)).*beta1.*lambda2)/(1+mo.*lambda2)+(lambda3 - lambda1).*((1-u1(j)).*beta2.*lambda3)/(1+m1.*lambda3)+lambda1.*(pi-mu-((1-u3(j)).*(beta2.*lambda4)/(lambda4+Ko)));
n12=G0+(lambda2 - lambda1).*(((1-u1(j)).*lambda1.*beta1)./(1+mo.*lambda2).^2)-lambda2.*(d1+(1-u2(j)).*eta1+(1-u3(j)).*alpha1)+lambda5.*(1-u2(j)).*eta1;
n13=G1+(lambda3-lambda1).*(((1-u1(j)).*lambda1.*omega.*beta2)./(1+m1.*lambda3).^2)-lambda3.*(d2+(1-u2(j)).*eta2+(1-u3(j)).*alpha2)+lambda5.*(1-u2(j)).*eta2;
n14=G2+(lambda2.*(1-(1-u3(j)).*beta3.*lambda1)./(lambda4+Ko))-lambda4.*delta1;
n15=-mu.*lambda5;
n21=((lambda2(j)-h2.*n12)-(lambda1(j)-h2.*n11)).*((1-u1(j)).*beta1.*lambda2(j))./(1+mo.*lambda2(j))+((lambda3(j)-h2.*n13)-(lambda1(j)-h2.*n11)).*((1-u1(j)).*beta2.*lambda3(j))/(1+m1.*lambda3(j))+(lambda1(j)-h2.*n11).*(pi-mu-((1-u3(j)).*(beta2.*lambda4(j))/(lambda4(j)+Ko)));
n22=G0+((lambda2(j)-h2.*n12)-(lambda1(j)-h2.*n11)).*(((1-u1(j)).*(lambda1(j)-h2.*n11).*beta1)./(1+mo.*x2).^2)-(lambda2(j)-h2.*n12).*(d1+(1-u2(j)).*eta1+(1-u3(j)).*alpha1)+(lambda5(j)-h2.*n15).*(1-u2(j)).*eta1;
n23=G1+((lambda3(j)-h2.*n13)-(lambda1(j)-h2.*n11)).*(((1-u1(j)).*(lambda1(j)-h2.*n11).*omega.*beta2)./(1+m1.*x3).^2)-(lambda3(j)-h2.*n13).*(d2+(1-u2(j)).*eta2+(1-u3(j)).*alpha2)+(lambda5(j)-h2.*n15).*(1-u2(j)).*eta2;
n24=G2+((lambda2(j)-h2.*n12).*(1-(1-u3(j)).*beta3.*(lambda1(j)-h2.*n11))./(lambda4(j)+Ko))-lambda4.*delta1;
n25=-mu.*(lambda5(j)-h2.*n15);
n31=((lambda2(j)-h2.*n22)-(lambda1(j)-h2.*n21)).*((1-u1(j)).*beta1.*x2)./(1+mo.*x2)+((lambda3(j)-h2.*n23)-(lambda1(j)-h2.*n21)).*((1-u1(j)).*beta2.*x3)/(1+m1.*x3)+(lambda1(j)-h2.*n21).*(pi-mu-((1-u3(j)).*(beta2.*x4)/(x4+Ko)));
n32=G0+((lambda2(j)-h2.*n22)-(lambda1(j)-h2.*n21)).*(((1-u1(j)).*(lambda1(j)-h2.*n21).*beta1)./(1+mo.*x2).^2)-(lambda2(j)-h2.*n22).*(d1+(1-u2(j)).*eta1+(1-u3(j)).*alpha1)+(lambda5(j)-h2.*n25).*(1-u2(j)).*eta1;
n33=G1+((lambda3(j)-h2.*n23)-(lambda1(j)-h2.*n21)).*(((1-u1(j)).*(lambda1(j)-h2.*n21).*omega.*beta2)./(1+m1.*x3).^2)-(lambda3(j)-h2.*n23).*(d2+(1-u2(j)).*eta2+(1-u3(j)).*alpha2)+(lambda5(j)-h2.*n25).*(1-u2(j)).*eta2;
n34=G2+((lambda2(j)-h2.*n22).*(1-(1-u3(j)).*beta3.*(lambda1(j)-h2.*n21))./(lambda4+Ko))-lambda4.*delta1;
n35=-mu.*(lambda5(j)-h2.*n25);
n41=((lambda2(j)-h2.*n32)-(lambda1(j)-h2.*n31)).*((1-u1(j)).*beta1.*x2)./(1+mo.*x2)+((lambda3(j)-h2.*n33)-(lambda1(j)-h2.*n31)).*((1-u1(j)).*beta2.*x3)/(1+m1.*x3)+(lambda1(j)-h2.*n31).*(pi-mu-((1-u3(j)).*(beta2.*x4)/(x4+Ko)));
n42=G0+((lambda2(j)-h2.*n32)-(lambda1(j)-h2.*n31)).*(((1-u1(j)).*(lambda1(j)-h2.*n31).*beta1)./(1+mo.*x2).^2)-(lambda2(j)-h2.*n32).*(d1+(1-u2(j)).*eta1+(1-u3(j)).*alpha1)+(lambda5(j)-h2.*n35).*(1-u2(j)).*eta1;
n43=G1+((lambda3(j)-h2.*n33)-(lambda1(j)-h2.*n31)).*(((1-u1(j)).*(lambda1(j)-h2.*n31).*omega.*beta2)./(1+m1.*x3).^2)-(lambda3(j)-h2.*n33).*(d2+(1-u2(j)).*eta2+(1-u3(j)).*alpha2)+(lambda5(j)-h2.*n35).*(1-u2(j)).*eta2;
n44=G2+((lambda2(j)-h2.*n32).*(1-(1-u3(j)).*beta3.*(lambda1(j)-h2.*n31))./(lambda4+Ko))-lambda4.*delta1;
n45=-mu.*(lambda5(j)-h2.*n35);
n51=((lambda2(j)-h2.*n42)-(lambda1(j)-h2.*n41)).*((1-u1(j)).*beta1.*x2)./(1+mo.*x2)+((lambda3(j)-h2.*n43)-(lambda1(j)-h2.*n41)).*((1-u1(j)).*beta2.*x3)/(1+m1.*x3)+(lambda1(j)-h2.*n41).*(pi-mu-((1-u3(j)).*(beta2.*x4)/(x4+Ko)));
n52=G0+((lambda2(j)-h2.*n42)-(lambda1(j)-h2.*n41)).*(((1-u1(j)).*(lambda1(j)-h2.*n41).*beta1)./(1+mo.*x2).^2)-(lambda2(j)-h2.*n42).*(d1+(1-u2(j)).*eta1+(1-u3(j)).*alpha1)+(lambda5(j)-h2.*n45).*(1-u2(j)).*eta1;
n53=G1+((lambda3(j)-h2.*n43)-(lambda1(j)-h2.*n41)).*(((1-u1(j)).*(lambda1(j)-h2.*n41).*omega.*beta2)./(1+m1.*x3).^2)-(lambda3(j)-h2.*n43).*(d2+(1-u2(j)).*eta2+(1-u3(j)).*alpha2)+(lambda5(j)-h2.*n45).*(1-u2(j)).*eta2;
n54=G2+((lambda2(j)-h2.*n42).*(1-(1-u3(j)).*beta3.*(lambda1(j)-h2.*n41))./(lambda4+Ko))-lambda4*delta1;
n55=-mu.*(lambda5(j)-h2.*n45);
dbstop if error
lambda1(j-1) = lambda1(j) - h./6.*(m11 + 2*m21 + 2*m31 + 2.*m41 + m51);
lambda2(j-1) = lambda2(j) - h/6*(m12 + 2*m22 + 2*m32 + 2.*m42+m52);
lambda3(j-1) = lambda3(j) - h/6*(m13 + 2*m23 + 2*m33 + 2.*m43+m53);
lambda4(j-1) = lambda4(j) - h/6*(m14 + 2*m24 + 2*m34 + 2.*m44+m54);
lambda5(j-1) = lambda5(j) - h/6.*(m15 + 2.*m25 + 2.*m35 + 2.*m55+m55);
u1 = min(0.9,max(0,((lambda2-lambda1).*(beta1.*lambda1.*lambda2)/(1+mo.*lambda2)+(lambda2-lambda1).*(beta1.*lambda3.*lambda1)/(1+m1.*lambda3))));
u2 = min(0.9,max(0,((lambda5-lambda2).*eta1.*lambda2+(lambda1-lambda3).*eta2.*lambda3)));
u3 = min(0.9,max(0,(lambda1-lambda2).*alpha1.*lambda2+(lambda4-lambda3).*alpha1.*lambda3-lambda1.*(beta3.*lambda1.*lambda4)/(lambda4+Ko)));
ua = 0.5.*(u1 + oldu1);
ub = 0.5.*(u2 + oldu2);
uc = 0.5.*(u3 + oldu3);
J=G0.*x2+G1.*x3+G2.*x5+0.5.*(K1.*u1.^2+K2.*u2.^2+K3.*u3.^2);
end
end
