Index in position 2 exceeds array bounds (must not exceed 1).

Hi. First of all sorry for asking the same question again and again but I am really close to the solution (probably).

Please do not mind long lines because most of them are constants and loops.

This code tries to solve 6 ODEs with 6 state variables [horizontal position (x1 and x2), altitude (x3), the true airspeed (x4), the heading angle (x5) and the mass of the aircraft (x6)] and 3 control inputs [engine thrust (u1), the bank angle (u2) and the flight path angle (u3)] by using Euler's Method.

Different flight maneuvers are performed for the specified time intervals.

Velocities.m, Cruise_Vel.m, Des_Vel.m, Thr_cl.m, Thr_cr.m, Thr_des.m, fuel_cl.m, fuel_cr.m, fuel_des.m,den.m,den_cr.m,den_des.m,drag.m,drag_cr.m,drag_des.m,lift.m,lift_cr.m,lift_des.m are functions in seperate tabs.

Main code is:

% Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
% Perform cruise flight for t=60 minutes.
% Turn with β=30 bank angle until heading is changed by η=270◦.
% Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
% Complete a 360◦ turn (loiter) at level flight.
% Descent to h3=800 [m] with κ=4.5◦ flight path angle.
% Aircraft Properties
W = .44225E+06;                             % .44225E+03 tons = .44225E+06 kg 
S = .51097E+03;                             % Surface Area [m^2]
g0 = 9.80665;                               % Gravitational acceleration [m/s2]

% solving 1st order ODE using numerical methods 
t0=0;           
tend=3960;
h=0.05; 
N=(tend-t0)/h;
t=t0:h:tend; 

% Preallocations
x = zeros(6,length(t));
x1 = zeros(1,length(t));
x2 = zeros(1,length(t));
x3 = zeros(1,length(t));
x4 = zeros(1,length(t));
x5 = zeros(1,length(t));
x6 = zeros(1,length(t));
u1 = zeros(1,length(t));
u2 = zeros(1,length(t));
u3 = zeros(1,length(t));
C_D= zeros(1,length(t));
p = zeros(1,length(t));
Cl = zeros(1,length(t));
f = zeros(1,length(t));
dx1dt = zeros(1,length(t));
dx2dt = zeros(1,length(t));
dx3dt = zeros(1,length(t));
dx4dt = zeros(1,length(t));
dx5dt = zeros(1,length(t));
dx6dt = zeros(1,length(t));

% Initial conditions
x(:,1)=[0;0;3608.92;1.0e+02 * 1.161544478045788;0;W];
for i=2:length(t)
      
if and (t(1,i-1) >= 0,t(1,i-1)<60)         % Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
   x3 = linspace(3608.92,5249.3,79201);
   x4 = Velocities(x3);                   % Changing speed [m/s]
   x5 = 0;                                           % Changing head angle [deg]
   f = fuel_cl(x3);                       % Changing fuel flow [kg/min] 
   
   u1 = Thr_cl(x3);                       % Changing thrust [N]
   u2 = 0;                                           % Changing bank angle [deg]
   u3 = 5;                                           % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag(x3,x4);                     % Changing drag coefficient
   Cl = lift(x3,x4);                      % Changing lift coefficient
   p = den(x3);                           % Changing density [kg/m3]                              
elseif and (t(1,i-1) >= 60,t(1,i-1)<3660) % Perform cruise flight for t=60 minutes.
   
   x3 = 5249.3;
   x4 = Cruise_Vel(x3);                          % Changing speed [m/s]
   x5 = 0;                                           % Changing head angle [deg]
   f = fuel_cr(x3);                              % Changing fuel flow [kg/min] 
   
   u1 = Thr_cr(x3);                              % Changing thrust [N]
   u2 = 0;                                           % Changing bank angle [deg]
   u3 = 0;                                           % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag_cr(x3,x4);                         % Changing drag coefficient
   Cl = lift_cr(x3,x4);                          % Changing lift coefficient
   p = den_cr(x3);                               % Changing density [kg/m3]
   
   elseif and (t(1,i-1) >= 3660,t(1,i-1)<3720) % Turn with β=30 bank angle until heading is changed by η=270◦.
   x3 = 5249.3;
   x4 = Cruise_Vel(x3);                          % Changing speed [m/s]
   x5 = 0:30:270;                                    % Changing head angle [deg]
   f = fuel_cr(x3);                              % Changing fuel flow [kg/min] 
   
   u1 = Thr_cr(x3);                              % Changing thrust [N]
   u2 = 30;                                          % Changing bank angle [deg]
   u3 = 0;                                           % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag_cr(x3,x4);                         % Changing drag coefficient
   Cl = lift_cr(x3,x4);                          % Changing lift coefficient
   p = den_cr(x3);                               % Changing density [kg/m3]
   
elseif and (t(1,i-1) >= 3720,t(1,i-1)<3780) % Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
   
   x3 = linspace(5249.3,3608.92,79201);
   x4 = Des_Vel(x3);                      % Changing speed [m/s]
   x5 = 270;                                         % Changing head angle [deg]
   f = fuel_des(x3);                      % Changing fuel flow [kg/min] 
   
   u1 = Thr_des(x3);                      % Changing thrust [N]
   u2 = 0;                                           % Changing bank angle [deg]
   u3 = 4;                                           % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag_des(x3,x4);                 % Changing drag coefficient
   Cl = lift_des(x3,x4);                  % Changing lift coefficient
   p = den_des(x3);                       % Changing density [kg/m3]
elseif and (t(1,i-1) >= 3780,t(1,i-1)<3900)      % Complete a 360◦ turn (loiter) at level flight.
   
   x3 = 3608.9;
   x4 = Cruise_Vel(x3);                          % Changing speed [m/s]
   lon = [270 300 360 60 120 180 240 270];
   x5 = wrapTo360(lon);                              % Changing head angle [deg]
   f = fuel_cr(x3);                              % Changing fuel flow [kg/min] 
   
   u1 = Thr_cr(x3);                              % Changing thrust [N] 
   u2 = 0;                                           % Changing bank angle [deg]
   u3 = 0;                                           % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag_cr(x3,x4);                         % Changing drag coefficient
   Cl = lift_cr(x3,x4);                          % Changing lift coefficient
   p = den_cr(x3);                               % Changing density [kg/m3]
elseif and (t(1,i-1) >= 3900,t(1,i-1)<3960)  % Descent to h3=800 [m] with κ=4.5◦ flight path angle.
  x3 = linspace(3608.92,2624.67,79201);
   
   x4 = Des_Vel(x3);                     % Changing speed [m/s] 
   x5 = 270;                                         % Changing head angle [deg]
   f = fuel_des(x3);                     % Changing fuel flow [kg/min] 
   
   u1 = Thr_des(x3);                     % Changing thrust [N] 
   u2 = 0;                                           % Changing bank angle [deg]
   u3 = 4.5;                                         % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag_des(x3,x4);                % Changing drag coefficient
   Cl = lift_des(x3,x4);                 % Changing lift coefficient
   p = den_des(x3);                      % Changing density [kg/m3]
else
    fprintf("A problem occured.");
end
dx1dt = x4 .* cos(x5) .* cos(u3);       
dx2dt = x4 .* sin(x5) .* cos(u3);       
dx3dt = x4 .* sin(u3);               
dx4dt = -C_D.*S.*p.*(x4.^2)./(2.*x6)-g0.*sin(u3)+u1./x6; 
dx5dt = -Cl.*S.*p.*x4./(2.*x6).*sin(u2);
dx6dt = -f;

x(1,i)= x(1,i-1) + h * dx1dt(1,i-1);     %%%%%%%%% line 138 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                 
x(2,i)= x(2,i-1) + h * dx2dt(1,i-1);
x(3,i)= x(3,i-1) + h * dx3dt(1,i-1);
x(4,i)= x(4,i-1) + h * dx4dt(1,i-1);
x(5,i)= x(5,i-1) + h * dx5dt(1,i-1);
x(6,i)= x(6,i-1) + h * dx6dt(1,i-1);
end
tot=cell2mat(f);                % Total fuel consumption during mission [kg/min]
Tot_fuel=sum(tot);
figure(1)
plot3(x1(:),x2(:),x3(:));       % 3D position graph
figure(2)
plot(t,x4(:));                  % Vtas − Time graph
figure(3)
plot(t,V_ver(:));               % V_vertical − Time graph
figure(4)
plot(t,x5(:));                  % Heading − Time graph
figure(5)
plot(t,x6(:));                  % Mass − Time graph
figure(6)
plot(t,u1(:));                  % Thrust − Time graph
figure(7)
plot(t,u2(:));                  % Bank Angle − Time graph
figure(8)
plot(t,u3(:));                  % Flight Path Angle − Time graph
fprintf('Total fuel consumption during mission is %.2f [kg]',Tot_fuel*tend/60);

The reason why I used 79201 sized array is length(t) = 79201.

And when I run:

Index in position 2 exceeds array bounds (must not exceed 1).

Error in forum (line 138)
x(1,i)= x(1,i-1) + h * dx1dt(1,i-1);

What should I do? Thanks.

One of the functions is seperate tabs is below, rest is similar:

function [Vtas_cl] = Velocities(x3)
%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Vcl_1 = 335;             % Standard calibrated airspeed [kt]
Vcl_2 = 172.3;           % Standard calibrated airspeed [kt] -> [m/s] (To find the  Mach transition altitude)
Vcl_2_in_knots = 335;    % Standard calibrated airspeed [kt] (To find the result in knots, if altitude is between 10,000 ft and Mach transition altitude)
M_cl  = 0.86;            % Standard calibrated airspeed [kt]

K = 1.4;                 % Adiabatic index of air
R = 287.05287;           % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065;           % ISA temperature gradient with altitude below the tropopause [K/m]
T0 = 288.15;             % Standard atmospheric temperature at MSL [K]
g0 = 9.80665;            % Gravitational acceleration [m/s2]
a0= 340.294;             % Speed of Sound [m/s]

Vd_CL1 = 5;              % Climb speed increment below 1500 ft (jet)
Vd_CL2 = 10;             % Climb speed increment below 3000 ft (jet)
Vd_CL3 = 30;             % Climb speed increment below 4000 ft (jet)
Vd_CL4 = 60;             % Climb speed increment below 5000 ft (jet)
Vd_CL5 = 80;             % Climb speed increment below 6000 ft (jet)

CV_min = 1.3;            % Minimum speed coefficient 
Vstall_TO = .14200E+03;  % Stall speed at take-off [KCAS]

CAS_climb  = Vcl_2;   
Mach_climb = M_cl;
delta_trans = (((1+((K-1)/2)*(CAS_climb/a0)^2)^(K/(K-1)))-1)/(((1+(K-1)/2*Mach_climb^2)^(K/(K-1)))-1);    % Pressure ratio at the transition altitude

teta_trans = delta_trans ^ (-Bt*R/g0);      % Temperature ratio at the transition altitude

H_p_trans_climb = (1000/0.348/6.5)*(T0*(1-teta_trans));   %  Transition altitude for climb [ft]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H_climb = x3; %%%%%% Input %%%%%%%%%%%%%%%%%%%
Vnom_climb_jet = zeros(1, length(H_climb));
for k1 = 1:length(H_climb)
    
    if (0<=H_climb(k1)&&H_climb(k1)<1500)
Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL1;

elseif (1500<=H_climb(k1)&&H_climb(k1)<3000)
        Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL2;
        
elseif (3000<=H_climb(k1)&&H_climb(k1)<4000)
        Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL3; 
        
elseif (4000<=H_climb(k1)&&H_climb(k1)<5000)
        Vnom_climb_jet (k1)= CV_min * Vstall_TO + Vd_CL4;
        
elseif (5000<=H_climb(k1)&&H_climb(k1)<6000)
        Vnom_climb_jet(k1) = CV_min * Vstall_TO + Vd_CL5;
        
elseif (6000<=H_climb(k1)&&H_climb(k1)<10000)
        Vnom_climb_jet (k1)= min(Vcl_1,250);
        
elseif (10000<=H_climb(k1)&&H_climb(k1)<=H_p_trans_climb)
        Vnom_climb_jet(k1) = Vcl_2_in_knots;
        
elseif (H_p_trans_climb<H_climb(k1))
        Vnom_climb_jet(k1) = M_cl;
    end

Vcas_cl(k1) = Vnom_climb_jet(k1)* 0.514;   % [kn] -> [m/s]


H_climb (k1)= H_climb(k1) * 0.3048;         % [feet] -> [m]
K = 1.4;                                    % Adiabatic index of air
R = 287.05287;                              % Real gas constant for air [m2/(K·s2)]
Bt = - 0.0065;                              % ISA temperature gradient with altitude below the tropopause [K/m]
deltaT = 0;                                 % Value of the real temperature T in ISA conditions [K]
T0 = 288.15;                                % Standard atmospheric temperature at MSL [K]
P0 = 101325;                                % Standard atmospheric pressure at MSL [Pa]
g0 = 9.80665;                               % Gravitational acceleration [m/s2]
p0 = 1.225;                                 % Standard atmospheric density at MSL [kg/m3]
visc = (K-1)./K;

T(k1) = T0 + deltaT + Bt * H_climb(k1);         % Temperature [K]
P (k1)= P0*((T(k1)-deltaT)/T0).^((-g0)/(Bt*R));     % Pressure [Pa]
p (k1)= P(k1) ./ (R*T(k1));                             % Density [kg/m^3]
Vtas_cl(k1) = (2*P(k1)/visc/p(k1)*((1 + P0/P(k1)*((1 + visc*p0*Vcas_cl(k1)*Vcas_cl(k1)/2/P0).^(1/visc)-1)).^(visc)-1)).^(1/2); % True Air Speed [m/s]
end
% Output
end

[Torsten]

   x4 = Velocities(x(3,i-1));
   x5 = 0;                                      % Changing head angle [deg]
   f = fuel_cl(x(3,i-1));                       % Changing fuel flow [kg/min] 
   
   u1 = Thr_cl(x(3,i-1));                       % Changing thrust [N]
   u2 = 0;                                      % Changing bank angle [deg]
   u3 = 5;                                      % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                          % Changing vertical speed [m/s]
   C_D = drag(x(3,i-1),x4);                     % Changing drag coefficient
   Cl = lift(x(3,i-1),x4);                      % Changing lift coefficient
   p = den(x(3,i-1));                           % Changing density [kg/m3]   

[Turgut Ataseven]

@Torsten Thanks for the response but I got the same error:

for i=2:length(t)
      
if and (t(1,i-1) >= 0,t(1,i-1)<60)              % Climb from h1=1100 [m] to h2=1600 [m] with α=5 flight path angle.
   x3 = linspace(3608.92,5249.3,79201);         % Changing altitude [m] -> [feet] 
   x4 = Velocities(x(3,i-1));                   % Changing speed [m/s]
   x5 = 0;                                      % Changing head angle [deg]
   f = fuel_cl(x(3,i-1));                       % Changing fuel flow [kg/min] 
   
   u1 = Thr_cl(x(3,i-1));                       % Changing thrust [N]
   u2 = 0;                                      % Changing bank angle [deg]
   u3 = 5;                                      % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                          % Changing vertical speed [m/s]
   C_D = drag(x(3,i-1),x4);                     % Changing drag coefficient
   Cl = lift(x(3,i-1),x4);                      % Changing lift coefficient
   p = den(x(3,i-1));                           % Changing density [kg/m3]                              
elseif and (t(1,i-1) >= 60,t(1,i-1)<3660) % Perform cruise flight for t=60 minutes.
   
   x3 = 5249.3;
   x4 = Cruise_Vel(x(3,i-1));                          % Changing speed [m/s]
   x5 = 0;                                           % Changing head angle [deg]
   f = fuel_cr(x(3,i-1));                              % Changing fuel flow [kg/min] 
   
   u1 = Thr_cr(x(3,i-1));                              % Changing thrust [N]
   u2 = 0;                                           % Changing bank angle [deg]
   u3 = 0;                                           % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag_cr(x(3,i-1),x4);                         % Changing drag coefficient
   Cl = lift_cr(x(3,i-1),x4);                          % Changing lift coefficient
   p = den_cr(x(3,i-1));                               % Changing density [kg/m3]
   
   elseif and (t(1,i-1) >= 3660,t(1,i-1)<3720) % Turn with β=30 bank angle until heading is changed by η=270◦.
   x3 = 5249.3;
   x4 = Cruise_Vel(x(3,i-1));                          % Changing speed [m/s]
   x5 = 0:30:270;                                    % Changing head angle [deg]
   f = fuel_cr(x(3,i-1));                              % Changing fuel flow [kg/min] 
   
   u1 = Thr_cr(x(3,i-1));                              % Changing thrust [N]
   u2 = 30;                                          % Changing bank angle [deg]
   u3 = 0;                                           % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag_cr(x(3,i-1),x4);                         % Changing drag coefficient
   Cl = lift_cr(x(3,i-1),x4);                          % Changing lift coefficient
   p = den_cr(x(3,i-1));                               % Changing density [kg/m3]
   
elseif and (t(1,i-1) >= 3720,t(1,i-1)<3780) % Descent from h2=1600 [m] to h1=1100 [m] with ζ=4◦ flight path angle.
   
   x3 = linspace(5249.3,3608.92,79201);
   x4 = Des_Vel(x(3,i-1));                      % Changing speed [m/s]
   x5 = 270;                                         % Changing head angle [deg]
   f = fuel_des(x(3,i-1));                      % Changing fuel flow [kg/min] 
   
   u1 = Thr_des(x(3,i-1));                      % Changing thrust [N]
   u2 = 0;                                           % Changing bank angle [deg]
   u3 = 4;                                           % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag_des(x(3,i-1),x4);                 % Changing drag coefficient
   Cl = lift_des(x(3,i-1),x4);                  % Changing lift coefficient
   p = den_des(x(3,i-1));                       % Changing density [kg/m3]
elseif and (t(1,i-1) >= 3780,t(1,i-1)<3900)      % Complete a 360◦ turn (loiter) at level flight.
   
   x3 = 3608.9;
   x4 = Cruise_Vel(x(3,i-1));                          % Changing speed [m/s]
   lon = [270 300 360 60 120 180 240 270];
   x5 = wrapTo360(lon);                              % Changing head angle [deg]
   f = fuel_cr(x(3,i-1));                              % Changing fuel flow [kg/min] 
   
   u1 = Thr_cr(x(3,i-1));                              % Changing thrust [N] 
   u2 = 0;                                           % Changing bank angle [deg]
   u3 = 0;                                           % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag_cr(x(3,i-1),x4);                         % Changing drag coefficient
   Cl = lift_cr(x(3,i-1),x4);                          % Changing lift coefficient
   p = den_cr(x(3,i-1));                               % Changing density [kg/m3]
elseif and (t(1,i-1) >= 3900,t(1,i-1)<3960)  % Descent to h3=800 [m] with κ=4.5◦ flight path angle.
  x3 = linspace(3608.92,2624.67,79201);
   
   x4 = Des_Vel(x(3,i-1));                     % Changing speed [m/s] 
   x5 = 270;                                         % Changing head angle [deg]
   f = fuel_des(x(3,i-1));                     % Changing fuel flow [kg/min] 
   
   u1 = Thr_des(x(3,i-1));                     % Changing thrust [N] 
   u2 = 0;                                           % Changing bank angle [deg]
   u3 = 4.5;                                         % Changing flight path angle [deg]
   V_ver = x4*sin(u3);                               % Changing vertical speed [m/s]
   C_D = drag_des(x(3,i-1),x4);                % Changing drag coefficient
   Cl = lift_des(x(3,i-1),x4);                 % Changing lift coefficient
   p = den_des(x(3,i-1));                      % Changing density [kg/m3]
else
    fprintf("A problem occured.");
end
dx1dt = x4 .* cos(x5) .* cos(u3);       
dx2dt = x4 .* sin(x5) .* cos(u3);       
dx3dt = x4 .* sin(u3);               
dx4dt = -C_D.*S.*p.*(x4.^2)./(2.*x6)-g0.*sin(u3)+u1./x6; 
dx5dt = -Cl.*S.*p.*x4./(2.*x6).*sin(u2);
dx6dt = -f;

x(1,i)= x(1,i-1) + h * dx1dt(1,i-1);            %%%%%%%% line 138 %%%%%%%%%%%%%%                                  
x(2,i)= x(2,i-1) + h * dx2dt(1,i-1);
x(3,i)= x(3,i-1) + h * dx3dt(1,i-1);
x(4,i)= x(4,i-1) + h * dx4dt(1,i-1);
x(5,i)= x(5,i-1) + h * dx5dt(1,i-1);
x(6,i)= x(6,i-1) + h * dx6dt(1,i-1);
end

After running:

Index in position 2 exceeds array bounds (must not exceed 1).

Error in Flight (line 138)
x(1,i)= x(1,i-1) + h * dx1dt(1,i-1);

I have tried to change x3 with this but it did not work either:

x(3,i-1) = linspace(3608.92,5249.3,79201);

Maybe x4 in C_D,Cl and p is the problematic one.

C_D = drag(x(3,i-1),x4);                     % Changing drag coefficient

[Torsten]

Up to which time t does the integration proceed now ?

[Turgut Ataseven]

@Torsten Hi, I get the error right after I run the code.

1x79201 double t is created. Some of x1,x2,...,dx1dt,dx2dt are 1x79201 doubles and some of them are just one single value.

[Turgut Ataseven]

i = 3 is interesting.

[Torsten]

I wonder where you define x6. I can't find it in the code you posted.

Further settings that will cause errors for later times:

x5 = 0:30:270;                                    % Changing head angle [deg]
lon = [270 300 360 60 120 180 240 270];

[Turgut Ataseven]

@Torsten Yes my bad. x6 is W and it is only defined at initial condition, there is not a formula for it but it won't stay constant.

Should have I gone with x(1,i) x(2,i)...x(6,i) dxdt(1,i)... dxdt(6,i) all the way along?

[Torsten]

I don't understand your question.

You don't need to save dx1dt,...,dx6dt. So no need to write them in arrays dx1dt(1,i-1),...,dx6dt(1,i-1). They can remain scalars.

Answer 1

[Jan]

dx1dt = x4 .* cos(x5) .* cos(u3);      % Now dx1dt is a scalar
...
x(1,i) = x(1,i-1) + h * dx1dt(1,i-1);  % This treats dx1dt as a vector

I guess, you want:

dx1dt(1, i) = x4 .* cos(x5) .* cos(u3);
%     ^^^^^

By the way:

if and (t(1,i-1) >= 0,t(1,i-1)<60)
    ...                              
elseif and (t(1,i-1) >= 60,t(1,i-1)<3660)    

can be simpilified. All t >= 0 and after the first test, t < 60 is excluded already:

if t(1, i-1) < 60
    ...                              
elseif t(1,i-1) < 3660    

There are more problems, e.g.:

plot3(x1(:),x2(:),x3(:));

x1 was pre-allocated, but now values have been applied. x3 is redefined multiple times.

[Turgut Ataseven]

Thanks for the response.

" x3 is redefined multiple times." has opened my eyes, I can say.