%This is the building's stiffness matrix
K=[6 -6 0; -6 18 -12; 0 -12 30];
K-K' % Be sure the stiffness matrix is symmetric
[Vk Dk]=eig(K) % Be sure the matrix is positive definite
% This is the mass matrix
M=[1.1 0 0; 0 1.5 0; 0 0 2]; 
% M*x''+K*x=f(t)
% Convert to state space s'=A*s
% where s(1)=x(1)  s(2)=x'(1) s(3)=x(2) etc.
A=[0 1 0 0 0 0; % Here A is the companion matrix
    -6/1.1 0 6/1.1 0 0 0;
    0 0 0 1 0 0;
    6/1.5 0 -18/1.5 0 12/1.5 0;
    0 0 0 0 0 1;
    0 0 12/2 0 -30/2 0];
[V D]=eig(A); % compute the eigenvalues D and eigenenvectors V
%ICs=[1 0 0 0 0 0]'; % try an initial condition (just the top floor)
ICs=V(:,1)+V(:,2); % try just the lowest freq mode
k=V\ICs;
s=@(t) k(1)*V(:,1)*exp(D(1,1)*t)+k(2)*V(:,2)*exp(D(2,2)*t)...
         +k(3)*V(:,3)*exp(D(3,3)*t)+k(4)*V(:,4)*exp(D(4,4)*t)...
         +k(5)*V(:,5)*exp(D(5,5)*t)+k(6)*V(:,6)*exp(D(6,6)*t)

count=1; % for frames of a movie
width=2; % width of the building in the plots
y=0:0.01:1; % this is for smoothly plotting the beams between floors
h01=y.^2.*(3-2*y);% this is for smoothly plotting the beams between floors
figure(1)
%step through time to make the movie of the building vibrating
for t=0:.1:6*pi
    sp=real(s(t));
    hold off
    plot([sp(1) sp(1)+width],[3 3],'r','LineWidth',3);
    hold on
    plot([sp(3) sp(3)+width],[2 2],'r','LineWidth',3);
    plot([sp(5) sp(5)+width],[1 1],'r','LineWidth',3);
    plot([0 0+width],[0 0],'r','LineWidth',3);
    plot(h01*sp(5),y,'r');plot(h01*sp(5)+width,y,'r');
    plot(h01*(sp(3)-sp(5))+sp(5),y+1,'r');plot(h01*(sp(3)-sp(5))+sp(5)+width,y+1,'r');
    plot(h01*(sp(1)-sp(3))+sp(3),y+2,'r');plot(h01*(sp(1)-sp(3))+sp(3)+width,y+2,'r');
    axis([-3 5 -1 4])
    %text(sp(1),3,'floor 1')
    Mv(count)=getframe;
    count=count+1;
end
movie(Mv);