Integrate function by for loop

clear
Cs1=4; Cs2=5.6; Cl1=2; Cl2=2.8; ns1=1.5; ns2=-1.5; nl1=1.5; nl2=-1.5;
rho=1.22; A=10; V=44.72; l=20; n=1;Ua=20;
yaw=0:0.5:90;
t=1:1:60;
u=[4.16886273000000 -5.46160601200000 -2.21536928300000 -0.107926889000000 4.17995704900000 2.70294649300000 -3.57644741500000 2.00451722300000 0.881703288000000 4.62541932200000 1.63637915200000 5.73084754400000 0.961007086000000 0.350184264000000 1.37768570200000 -7.64386640500000 1.09384013000000 -3.32210770200000 1.06977546700000 3.77699477800000 -2.80761673400000 -2.51814521100000 -2.56997163700000 -0.434123529000000 -7.33456835300000 3.18164556400000 -2.46259887500000 -2.81705263700000 -6.11824449800000 4.20959983600000 0.958003023000000 1.34911158700000 -3.66536927800000 1.33702769200000 0.951579011000000 4.81211100200000 1.70504507900000 7.08168644200000 3.20839901100000 -0.438521707000000 2.92055007900000 -0.280320940000000 -1.25512163200000 -2.68273256600000 -5.07343734300000 -1.42260189600000 4.56907055100000 3.27623428200000 -1.23514645600000 -2.86948142600000 1.75847520500000 2.79193783800000 -4.18180552400000 -3.14073219500000 -2.43496557500000 -1.28117719800000 4.95079755900000 -3.52809836900000 1.96375591300000 -5.66592848200000]

for i=1:length(yaw)
ms(i)=2*sin(yaw(i)*pi/180)*V/l;
ml(i)=2.5*sin(yaw(i)*pi/180)*V/l;

if yaw(i)<=40
Cs(i)=Cs1*(sin(yaw(i)*pi/180)/sin(30*pi/180))^ns1;
Cl(i)=Cl1*(sin(yaw(i)*pi/180)/sin(30*pi/180))^nl1;
elseif yaw(i)>40
Cs(i)=Cs2*(sin(yaw(i)*pi/180)/sin(30*pi/180))^ns2;
Cl(i)=Cl2*(sin(yaw(i)*pi/180)/sin(30*pi/180))^nl2;
end

S1(i)=0.5*rho*A*(V^2)*Cs(i); % Steady side forces (Uniform along each row)
L1(i)=0.5*rho*A*(V^2)*Cl(i); % Steady lift forces (Uniform along each row)

Msi=ms(i);
Mli=ml(i);
Csi=Cs(i);
Cli=Cl(i);
Yawi=yaw(i*pi/180);

for j=1:length(t)
Tj=t(j);
Uj=u(j);
fun1=@(t1)(((2.*pi.*n.*Msi)^2).*t1.*exp(-2.*pi.*n.*Msi.*t1)).*Uj.*(Tj-t1);
fun2=@(t1)(((2.*pi.*n.*Mli)^2).*t1.*exp(-2.*pi.*n.*Mli.*t1)).*Uj.*(Tj-t1);
S2(i,j)=rho*A*Csi*Ua*(1+(1/(2*Csi))*diff(Csi)*cot(Yawi))*integral(fun1,0,inf); % Unsteady side forces (noniform along each row and column)
L2(i,j)=rho*A*Csi*Ua*(1+(1/(2*Csi))*diff(Csi)*cot(Yawi))*integral(fun2,0,inf); % Unsteady lift forces (noniform along each row and column)
S(i,j)=Sl(i)+S2(i,j); % Total side forces
L(i,j)=S1(i)+S2(i,j); % Total lift forces
j=j+1
end

i=i+1;
end
Dear experts:
The equations I am going to solve are attached as above screenshot and the code is listed.
To summary, the equation is the: Total force(i,j)=Mean force (i)+Unsteady force (i,j (involve intergration over time interval t1)). Hence I will get i*j number of total force result. It is 181*60 (row*column) data.
The mean force componant S1, L1 is uniform value inside each row, which mean that the mean force does not vary along columnn on each row (i); but the mean force componants only has different value between rows.
The calculation of the unsteady force S2,L2 is more complex and the unsteady force is different between each row and column (i,j), and inside each (i,j), the unsteady force is calculated based on the intergration over time interal (0,infinite):S2 (i,j)=(eq…)*integration(fun1(t1), 0 ,inf), L2 (i,j)=(eq…)*integration(fun2(t1), 0 ,inf),
Thereforce, the total force S(i,j)=S1(i)+S2(i,j), L(i,j)=L1(i)+L2(i,j)
The code does not have problem with calculating the mean force (i) S1, L1, but it shows "Array indices must be positive integers or logical values" after I add the for loop (j) to calculate the unsteay force (i,j) S2, L2 which involve integration over time interval t1.
Could our expert help me fix the code?
Kind Regardsclear
Cs1=4; Cs2=5.6; Cl1=2; Cl2=2.8; ns1=1.5; ns2=-1.5; nl1=1.5; nl2=-1.5;
rho=1.22; A=10; V=44.72; l=20; n=1;Ua=20;
yaw=0:0.5:90;
t=1:1:60;
u=[4.16886273000000 -5.46160601200000 -2.21536928300000 -0.107926889000000 4.17995704900000 2.70294649300000 -3.57644741500000 2.00451722300000 0.881703288000000 4.62541932200000 1.63637915200000 5.73084754400000 0.961007086000000 0.350184264000000 1.37768570200000 -7.64386640500000 1.09384013000000 -3.32210770200000 1.06977546700000 3.77699477800000 -2.80761673400000 -2.51814521100000 -2.56997163700000 -0.434123529000000 -7.33456835300000 3.18164556400000 -2.46259887500000 -2.81705263700000 -6.11824449800000 4.20959983600000 0.958003023000000 1.34911158700000 -3.66536927800000 1.33702769200000 0.951579011000000 4.81211100200000 1.70504507900000 7.08168644200000 3.20839901100000 -0.438521707000000 2.92055007900000 -0.280320940000000 -1.25512163200000 -2.68273256600000 -5.07343734300000 -1.42260189600000 4.56907055100000 3.27623428200000 -1.23514645600000 -2.86948142600000 1.75847520500000 2.79193783800000 -4.18180552400000 -3.14073219500000 -2.43496557500000 -1.28117719800000 4.95079755900000 -3.52809836900000 1.96375591300000 -5.66592848200000]

for i=1:length(yaw)
ms(i)=2*sin(yaw(i)*pi/180)*V/l;
ml(i)=2.5*sin(yaw(i)*pi/180)*V/l;

S1(i)=0.5*rho*A*(V^2)*Cs(i); % Steady side forces (Uniform along each row)
L1(i)=0.5*rho*A*(V^2)*Cl(i); % Steady lift forces (Uniform along each row)

Msi=ms(i);
Mli=ml(i);
Csi=Cs(i);
Cli=Cl(i);
Yawi=yaw(i*pi/180);

i=i+1;
end
Dear experts:
The equations I am going to solve are attached as above screenshot and the code is listed.
To summary, the equation is the: Total force(i,j)=Mean force (i)+Unsteady force (i,j (involve intergration over time interval t1)). Hence I will get i*j number of total force result. It is 181*60 (row*column) data.
The mean force componant S1, L1 is uniform value inside each row, which mean that the mean force does not vary along columnn on each row (i); but the mean force componants only has different value between rows.
The calculation of the unsteady force S2,L2 is more complex and the unsteady force is different between each row and column (i,j), and inside each (i,j), the unsteady force is calculated based on the intergration over time interal (0,infinite):S2 (i,j)=(eq…)*integration(fun1(t1), 0 ,inf), L2 (i,j)=(eq…)*integration(fun2(t1), 0 ,inf),
Thereforce, the total force S(i,j)=S1(i)+S2(i,j), L(i,j)=L1(i)+L2(i,j)
The code does not have problem with calculating the mean force (i) S1, L1, but it shows "Array indices must be positive integers or logical values" after I add the for loop (j) to calculate the unsteay force (i,j) S2, L2 which involve integration over time interval t1.
Could our expert help me fix the code?
Kind Regards clear
Cs1=4; Cs2=5.6; Cl1=2; Cl2=2.8; ns1=1.5; ns2=-1.5; nl1=1.5; nl2=-1.5;
rho=1.22; A=10; V=44.72; l=20; n=1;Ua=20;
yaw=0:0.5:90;
t=1:1:60;
u=[4.16886273000000 -5.46160601200000 -2.21536928300000 -0.107926889000000 4.17995704900000 2.70294649300000 -3.57644741500000 2.00451722300000 0.881703288000000 4.62541932200000 1.63637915200000 5.73084754400000 0.961007086000000 0.350184264000000 1.37768570200000 -7.64386640500000 1.09384013000000 -3.32210770200000 1.06977546700000 3.77699477800000 -2.80761673400000 -2.51814521100000 -2.56997163700000 -0.434123529000000 -7.33456835300000 3.18164556400000 -2.46259887500000 -2.81705263700000 -6.11824449800000 4.20959983600000 0.958003023000000 1.34911158700000 -3.66536927800000 1.33702769200000 0.951579011000000 4.81211100200000 1.70504507900000 7.08168644200000 3.20839901100000 -0.438521707000000 2.92055007900000 -0.280320940000000 -1.25512163200000 -2.68273256600000 -5.07343734300000 -1.42260189600000 4.56907055100000 3.27623428200000 -1.23514645600000 -2.86948142600000 1.75847520500000 2.79193783800000 -4.18180552400000 -3.14073219500000 -2.43496557500000 -1.28117719800000 4.95079755900000 -3.52809836900000 1.96375591300000 -5.66592848200000]

for i=1:length(yaw)
ms(i)=2*sin(yaw(i)*pi/180)*V/l;
ml(i)=2.5*sin(yaw(i)*pi/180)*V/l;

S1(i)=0.5*rho*A*(V^2)*Cs(i); % Steady side forces (Uniform along each row)
L1(i)=0.5*rho*A*(V^2)*Cl(i); % Steady lift forces (Uniform along each row)

Msi=ms(i);
Mli=ml(i);
Csi=Cs(i);
Cli=Cl(i);
Yawi=yaw(i*pi/180);

Cart

Cart