I keep receiving an error when I try to run this code for a 3DOF mass-spring system "Index in position 1 exceeds array bounds. Index must not exceed 1".I am confused as to why

I am trying to run this code and everytime I do, I always get an error message that says thus: "Index in position 1 exceeds array bounds. Index must not exceed 1". I am clearly doing something wrong and cannot seem to figure out what it is.

Here is my code below:

clear all;
close all;
n=3; %degree of freedom
m1= 45400; %mass for first floor in (kg)
m2= 45400; %mass for second floor in (kg)
m3= 45400; %mass for third floor in (kg)
k1= 120000; %first steel stiffness in (N/m)
k2= 120000; %second steel stiffness in (N/m)
k3= 120000; %third steel stiffness in (N/m)
c1= 7381.06; %first damper in (N.s/m^2)
c2= 7381.06; %second damper in (N.s/m^2)
c3= 7381.06; %third damper in (N.s/m^2)
M= [m1 0 0;0 m2 0;0 0 m3]; %mass matrix
K= [k1+k2,-k2,0;-k2,k2+k3,-k3;0,-k3,k3]; %stiffness matrix
C= [c1+c2,-c2,0;-c2,c2+c3,-c3;0,-c3,c3]; % damping matrix
[u,s]=eig(K,M);
u1=u;
s1=s;
x0=[0 0 0]';
xd0=[0 0 0]';
Kh=M^(-0.5)*K*M^(-0.5);
[P,s2]=eig(Kh);
w=sqrt(diag(s2));
S=M^(-0.5)*P;
r0=inv(S)*xd0;
rd0=inv(S)*x0;
%alp=0;
%bet=0.1;
zeta=0.05;
wd=w.*((1-zeta^2)).^0.5;
F=[0 0 840]';
f=P'*M^(-.5)*F;
wr=[0 0 0]';
XX=(f./(((w.^2-wr.^2)+((2*zeta.*w.*wr).^2)).^0.5));
th=atan2((2*zeta.*w.*wr),((w.^2)-(wr.^2)));
phir = atan2(((wd.*(r0-XX.*cos(th)))),(rd0+(zeta.*w.*(r0-XX.*cos(th))))-(wr.*XX.*sin(th)));
Ar=(r0-(XX.*cos(th)))./sin(phir);
t=(0:0.001:10)';
for i=1:n
    rt(i,:)=Ar(i,1).*exp(-zeta(i,1)*w(i,1)*t).*sin(wd(i,1)*t+phir(i,1))+XX(i,1)*cos(wr(i,1)*t-th(i,1));
end

Index in position 1 exceeds array bounds. Index must not exceed 1.

%xt=S*rt;
%xt=xt';
figure(1)
subplot(3,1,1)
plot(t,rt(1,:))
title ('Response of Mass #1');
xlabel('Time: in (seconds)')
ylabel('Response: in (meters)')
grid on
subplot(3,1,2)
plot(t,rt(1,:))
title ('Response of Mass #2');
xlabel('Time: in (seconds)')
ylabel('Response: in (meters)')
grid on
subplot(3,1,3)
plot(t,rt(1,:))
title ('Response of Mass #3');
xlabel('Time: in (seconds)')
ylabel('Response: in (meters)')
grid on

Answer 1

[Cris LaPierre]

The problem is that zeta is a scalar, not a vector. The fix is to not index this variable.

n=3; %degree of freedom
m1= 45400; %mass for first floor in (kg)
m2= 45400; %mass for second floor in (kg)
m3= 45400; %mass for third floor in (kg)
k1= 120000; %first steel stiffness in (N/m)
k2= 120000; %second steel stiffness in (N/m)
k3= 120000; %third steel stiffness in (N/m)
c1= 7381.06; %first damper in (N.s/m^2)
c2= 7381.06; %second damper in (N.s/m^2)
c3= 7381.06; %third damper in (N.s/m^2)
M= [m1 0 0;0 m2 0;0 0 m3]; %mass matrix
K= [k1+k2,-k2,0;-k2,k2+k3,-k3;0,-k3,k3]; %stiffness matrix
C= [c1+c2,-c2,0;-c2,c2+c3,-c3;0,-c3,c3]; % damping matrix
[u,s]=eig(K,M);
u1=u;
s1=s;
x0=[0 0 0]';
xd0=[0 0 0]';
Kh=M^(-0.5)*K*M^(-0.5);
[P,s2]=eig(Kh);
w=sqrt(diag(s2));
S=M^(-0.5)*P;
r0=inv(S)*xd0;
rd0=inv(S)*x0;
%alp=0;
%bet=0.1;
zeta=0.05;
wd=w.*((1-zeta^2)).^0.5;
F=[0 0 840]';
f=P'*M^(-.5)*F;
wr=[0 0 0]';
XX=(f./(((w.^2-wr.^2)+((2*zeta.*w.*wr).^2)).^0.5));
th=atan2((2*zeta.*w.*wr),((w.^2)-(wr.^2)));
phir = atan2(((wd.*(r0-XX.*cos(th)))),(rd0+(zeta.*w.*(r0-XX.*cos(th))))-(wr.*XX.*sin(th)));
Ar=(r0-(XX.*cos(th)))./sin(phir);
t=(0:0.001:10)';
for i=1:n
    rt(i,:)=Ar(i,1).*exp(-zeta*w(i,1)*t).*sin(wd(i,1)*t+phir(i,1))+XX(i,1)*cos(wr(i,1)*t-th(i,1));
    %                        ^^ removed (i,1)
end
%xt=S*rt;
%xt=xt';
figure(1)
subplot(3,1,1)
plot(t,rt(1,:))
title ('Response of Mass #1');
xlabel('Time: in (seconds)')
ylabel('Response: in (meters)')
grid on
subplot(3,1,2)
plot(t,rt(1,:))
title ('Response of Mass #2');
xlabel('Time: in (seconds)')
ylabel('Response: in (meters)')
grid on
subplot(3,1,3)
plot(t,rt(1,:))
title ('Response of Mass #3');
xlabel('Time: in (seconds)')
ylabel('Response: in (meters)')
grid on

[Fiyinfoluwa Adeseha]

Thank you so very much.