How do I rectify this 'out of index' error?

function [ s ] = harmonic(ms,omegas,damprat,ug,omegax)
ks=ms*omegas^2;
u(1)=0;
u(2)=0;
a(1)=0;
cd=2*ms*omegas*damprat;
k=ms/(0.01)^2+cd/0.02;
x=ms/(0.01)^2-cd/0.02;
y=ks-2*ms/(0.01)^2;
for t=0:0.01:10
    for n=1:1:1001
    if n==(t*100+1)
    p(n)=-ms*ug*sin(omegax*t)-x*u(n)-y*u(n+1);
    u(n+2)=p(n)/k;
    a(n+1)=(u(n+2)-2*u(n+1)+u(n))/(0.01)^2;
    fprintf('t=%5d\t%10.3f\n',t,a(n))
    plot(t,a(n))
    end
    end
end

When I run the above function in command window, I get the following:

harmonic(820.116,7.068,0.014,70,6.283)
t=    0              0.000
t=1.000000e-002      0.000
t=2.000000e-002     -4.352
t=3.000000e-002     -6.449
t=4.000000e-002     -5.280
t=5.000000e-002     -1.471
t=6.000000e-002      3.059
t=7.000000e-002      6.072
t=8.000000e-002      6.123
t=9.000000e-002      3.260
t=1.000000e-001     -1.027
t=1.100000e-001     -4.576
t=1.200000e-001     -5.631
t=1.300000e-001     -3.705
??? Attempted to access u(17); index out of bounds because numel(u)=16.

Error in ==> harmonic at 13
    p(n)=-ms*ug*sin(omegax*t)-x*u(n)-y*u(n+1);

Please suggest why does this happen?

Answer 1

[Walter Roberson]

http://matlab.wikia.com/wiki/FAQ#Why_is_0.3_-_0.2_-_0.1_.28or_similar.29_not_equal_to_zero.3F

[Maneet]

Thanks Mr. Walter for replying. Could you suggest anything for rectifying this error?

Answer 2

[Amit]

The error is as Walter mentioned. If you want to stick to the format you showed in your post, the solution is also included in the link posted by Walter (istead of using a == b, use (a-b) < tol or use n == ceil(t*100+1)

However, I am not sure why you are using so many loops to calculate what you are calculating? why dont you use something like this (this is not the most optimized code, but without knowing the context I dont wanna change it from original too much):

n = 1:1:101;
t = (n-1)/100;
u = zeros(103,1); % Initialize u
a = zeros(102,1);  % initialize a
for ii = 1:numel(n)
    p(ii)=-ms*ug*sin(omegax*t(ii))-x*u(ii)-y*u(ii+1);
    u(ii+2)=p(ii)/k;
    a(ii+1)=(u(ii+2)-2*u(ii+1)+u(ii))/(0.01)^2;
    fprintf('t=%5d\t%10.3f\n',t(ii),a(ii))
    plot(t(ii),a(ii))
end

Also, without hold the plot will not give anything. what you are trying to plot? A line or just data points?

[Maneet]

Thnks a lot, Amit. Using abs(a-b)<tol, worked. And will try your code as well. Thnks for your time.