clc
clear

k=1000
K=[2*k -k 0
    -k 2*k -k
    0 -k 2*k]
M=[10 0 0
    0 20 0
    0 0 10]

pause

[V,D]=eig(K,M)

mr=V'*M*V
kr=V'*K*V

omega=diag(sqrt(D))

pause

i=1;
j=3;

s=0;
 for W=0:.1:30
     HH=0;
     s=s+1;
for r=1:3
    HH=HH+V(r,i)*V(r,j)/(omega(r)^2-W^2);
end
H(s)=HH;
 end
W1=[0:.1:30];
plot(W1,20*log10(abs(H)),'r*')
pause
hold on



W1=[0:.1:30];
H1=zeros(size(W1));
for r=1:3
    H1=H1+V(r,i)*V(r,j)./(omega(r)^2-W1.^2);
end

plot(W1,20*log10(abs(H1)))