Loop not returning correct values

%% Clear MATLAB
clear;
clc;
%% Parameters
tlength=16.74; % time (s)
h=0.02; % time step (s)
N=1:ceil(tlength/h);
t=0.02:h:tlength; % time steps
w0=31.416; % Natural circular freq
zeta0=0.05;
eta=2;
tau=0.12;
m=7; 
R=eta*m*w0^2;
%% Initial Conditions
t(1)=0.02;
v(1)=0; % Velocity at t=0.02s
y3(1)=0; % Displacement of EDD at t=0s
ag(1)=-0.1; % Ground Acceleration at t=0s
a(1)=0; % Acceleration of frame
u(1)=1; % Displacement of frame
vedd(1)=0; % Velocity of EDD.
%% Define Equations
f1=@(t,u,v,y3) v;
f2=@(t,u,v,y3) -u*w0^2-2*zeta0*w0*v-eta*y3*w0^2-ag;
f3=@(t,u,v,y3) v-y3/tau;
%% RK4 Loop
for i=1:N
    t(i+1)=t(i)+h;
    
    F1_1=h*f1(t(i),u(i),v(i),y3(i));
    F2_1=h*f2(t(i),u(i),v(i),y3(i));
    F3_1=h*f3(t(i),u(i),v(i),y3(i));
    F1_2=h*f1(t(i)+h/2,u(i)+h/2*F1_1,v(i)+h/2*F1_1,y3(i)+h/2*F1_1);
    F2_2=h*f2(t(i)+h/2,u(i)+h/2*F2_1,v(i)+h/2*F2_1,y3(i)+h/2*F2_1);
    F3_2=h*f3(t(i)+h/2,u(i)+h/2*F3_1,v(i)+h/2*F3_1,y3(i)+h/2*F3_1);
    F1_3=h*f1(t(i)+h/2,u(i)+h/2*F1_2,v(i)+h/2*F1_2,y3(i)+h/2*F1_2);
    F2_3=h*f2(t(i)+h/2,u(i)+h/2*F2_2,v(i)+h/2*F2_2,y3(i)+h/2*F2_2);
    F3_3=h*f3(t(i)+h/2,u(i)+h/2*F3_2,v(i)+h/2*F3_2,y3(i)+h/2*F3_2);
    F1_4=h*f1(t(i)+h,u(i)+h*F1_3,v(i)+h*F1_3,y3(i)+h*F1_3);
    F2_4=h*f2(t(i)+h,u(i)+h*F2_3,v(i)+h*F2_3,y3(i)+h*F2_3);
    F3_4=h*f3(t(i)+h,u(i)+h*F3_3,v(i)+h*F3_3,y3(i)+h*F3_3);
    v(i+1)= v(i)+1/6*(F1_1+2*F1_2+2*F1_3+F1_4);
    a(i+1)= a(i)+1/6*(F2_1+2*F2_2+2*F2_3+F2_4);
    vedd(i+1)= vedd(i)+1/6*(F3_1+2*F3_2+2*F3_3+F3_4);
end

All that is being returned is v=[0,0], vedd=[0,0], and a=[0,-11.2904] but is should be returning all N values.

Answer 1

[Image Analyst]

For some weird reason you made N a vector instead of a scalar so the for loop line does not make sense. Make it a scalar:

N = ceil(tlength/h);   % NOT 1 : ceil(tlength/h) !!!
for i = 1 : N
    fprintf('i = %d of %d.\n', i, N)

[joe brady]

Ok - lower the tone.

Thanks for the reply all sorted now 👍