Index exceeds the number of array elements

When running the below code, my counter ii stops at 501 rather than stepping to 1552 like u. Tried the same approach as how u was loaded with Nan's but just run into different errors.

m = 290/386.4;                  % mass [kip s2/in]
L = 14*12;              %Height [in.]
E = 29000;              %Modulus of Elasticity [psi]
I = 209;                %Moment of Inertia [in4]
n = 24;                 %Number of Columns
Is = I*n;               %Total Moment of Inertia [in4]
k = 3*E*Is/(L^3);       %Stiffness [lb/in]
omegan = (k/m)^.5;      % natural frequency [rad/s]
Tn = 2*pi/omegan       % natural period [s]
zeta = 0.05;            % damping ratio []
omegad = omegan*sqrt(1-zeta^2); % damped natural frequency [rad/s]
c = 2*zeta*omegan*m;    % damping [lb*s/in]
u0 = 0;                 % initial displacement [in]
ud0 = 0;                % initial velocity [in/s]

Define the loading time history .

% Load ground motion -- El Centro
ugdd = load('ElCentro.dat');
ugdd = reshape(ugdd',[],1);             % ground acceleration [g]
Dt = 0.02;                              % time increment [s]
t = (0:length(ugdd)-1)*Dt;              % time [sec]
[ugddo,i_ugddo] = max(abs(ugdd));
figure(1);
subplot(411)
    plot(t,ugdd,t(i_ugddo),ugdd(i_ugddo),'o','linewidth',2)
    ylim([-.4,.4]); xlim([0,30])
    ylabel('u_g'''' [g]')
    title(['u"_{go} = ',num2str(max(abs(ugdd))),' g'])
    grid on

Central Difference Method

Initialize calculations for central difference method:

• Acceleration at t = 0:
• Displacement at :
• Coefficient on :
• Coefficient on :
• Coefficient on :

u = nan(size(t));               % allocate memory for u time history 
u(1) = u0;                      % displacement at t=0 [m]
p0 = p(1);                      % load at t=0 [N]
udd0 = (p0-c*ud0-k*u0)/m;       % acceleration at t=0 [m/s/s]
u_1 = u0 - Dt*ud0 + Dt^2/2*udd0;% displacement at t=-dt [m]
khat = m/Dt^2 + c/(2*Dt);       % coefficient on u_(i+1)
a  = m/Dt^2 - c/(2*Dt);         % coefficient on u_(i-1) 
b  = k - 2*m/Dt^2;              % coefficient on u_i

Iterate centeral difference method. For time step i:

1. 
2. 
3. If required: and

ii = 1;                                % initialize counter
u(ii+1) = (p(ii)-a*u_1-b*u(ii))/khat;  % first iteration
for ii = 2:length(t)-1                % loop over remaining time
    phat = p(ii)-a*u(ii-1)-b*u(ii);    % phat at i
    u(ii+1) = phat/khat;               % displacement at i+1
end

[Walter Roberson]

p0 = p(1);                      % load at t=0 [N]

p is not defined before that point, so we have no idea where it is coming from or how large it is. As outside observers we have no reason to expect that it is at least as long as u is.

[Daniel Murphy]

Understood, thank you