Error of Attempted to access u1(51); index out of bounds because numel(u1)=1. please the error is in line 111

Error of Attempted to access u1(51);

index out of bounds because numel(u1)=1.

Please help. The error is in line 111.

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

Answer 1

[Image Analyst]

Your "ul" did have 51 elements until you did this to set it equal to the min, which is a scalar:

u1 = min(0.9,max(0,((lambda2-lambda1).*(beta1.*lambda1.*lambda2)/(1+mo.*lambda2)+(lambda2-lambda1).*(beta1.*lambda3.*lambda1)/(1+m1.*lambda3))));

I have no idea what you're doing there but you're not getting a vector out of that so you should take a look at that again.

The very disappointing lack of comments prevents us from understanding what you're doing. Sorry but it just looks like alphabet soup to me so I'll have to leave it to you to figure out how to not turn ul from a vector into a scalar. I formatted your code for you to make it easy for people to copy and paste into MATLAB but that's about all the time I can spend on it. Maybe someone else would love to take a deeper dive into it than me so be very patient as you debug in parallel. If you solve it, let us know so nobody wastes time solving an already-solved problem.