close all; m=[1 1 1]'; k=[3 2 1]';       % masses and spring constants
                                         % rest positions and initial
l0=[1 1 1]'; u0=[0.1, 0.1, 0.9]';        % displacements
A=[(-k(1)-k(2))/m(1), k(2)/m(1),          0
   k(2)/m(2),         (-k(2)-k(3))/m(2), k(3)/m(2)
   0,                 k(3)/m(3),         -k(3)/m(3)];
[U,l]=eig(A); l=diag(l); sl=sqrt(l); Ux=[U,u0];
w=1;                                     % w is the box width
for sp=1:4                               % draw four plots:
  subplot(4,1,sp);                       % 1-3: eigensolutions
  axis([-1,10,0,3]);                     % 4: composite solution
  axis equal; axis off; hold on;
  o(1)=l0(1);                            % offsets for drawing
  o(2)=l0(1)+Ux(1,sp)+w+l0(2);
  o(3)=l0(1)+Ux(1,sp)+w+l0(2)+Ux(2,sp)+w+l0(3);
  bs(:,sp)=[                             % box plots
  plot([o(1)+Ux(1,sp),o(1)+Ux(1,sp)+w,o(1)+Ux(1,sp)+w,o(1)+ ...
    Ux(1,sp),o(1)+Ux(1,sp)],[0,0,w,w,0],'k-');    
  plot([o(2)+Ux(2,sp),o(2)+Ux(2,sp)+w,o(2)+Ux(2,sp)+w,o(2)+ ...
    Ux(2,sp),o(2)+Ux(2,sp)],[0,0,w,w,0],'k-');
  plot([o(3)+Ux(3,sp),o(3)+Ux(3,sp)+w,o(3)+Ux(3,sp)+w,o(3)+ ...
    Ux(3,sp),o(3)+Ux(3,sp)],[0,0,w,w,0],'k-'); ];
  ss(:,sp)=[                             % spring plots
  plot([0:(o(1)+Ux(1,sp))/10:o(1)+Ux(1,sp)], ... 
    [w/2+w/10,w/2+w/10*cos((1:9)*pi),w/2+w/10],'b-');
  plot([o(1)+Ux(1,sp)+w:(o(2)-o(1)+Ux(2,sp)-Ux(1,sp)-w)/10: ...
    o(2)+Ux(2,sp)],[w/2+w/10,w/2+w/10*cos((1:9)*pi),w/2+w/10],'b-');
  plot([o(2)+Ux(2,sp)+w:(o(3)-o(2)+Ux(3,sp)-Ux(2,sp)-w)/10: ...
    o(3)+Ux(3,sp)],[w/2+w/10,w/2+w/10*cos((1:9)*pi),w/2+w/10],'b-');
  ];
  tmp=plot([0,0,10],[2*w,0,0],'k-');     % draw axes
  set(tmp,'LineWidth',3);  
end
drawnow;                                 % draw initial positions
t=0; dt=0.04;                            % perform simulation in time
while 1 
  t=t+dt;                                % compute the new position
  v=real(U*diag(exp(sl*t))*(U\Ux));
  for sp=1:4
    o(1)=l0(1);                          % update the drawing
    o(2)=l0(1)+v(1,sp)+w+l0(2);
    o(3)=l0(1)+v(1,sp)+w+l0(2)+v(2,sp)+w+l0(3);
    set(bs(1,sp),'XData',[o(1)+v(1,sp),o(1)+v(1,sp)+ ...        
      w,o(1)+v(1,sp)+w,o(1)+v(1,sp),o(1)+v(1,sp)]);
    set(bs(2,sp),'XData',[o(2)+v(2,sp),o(2)+v(2,sp)+ ... 
      w,o(2)+v(2,sp)+w,o(2)+v(2,sp),o(2)+v(2,sp)]);
    set(bs(3,sp),'XData',[o(3)+v(3,sp),o(3)+v(3,sp)+ ... 
      w,o(3)+v(3,sp)+w,o(3)+v(3,sp),o(3)+v(3,sp)]);
    set(ss(1,sp),'XData',[0:(o(1)+v(1,sp))/10:o(1)+v(1,sp)]);
    set(ss(2,sp),'XData',[o(1)+v(1,sp)+w:(o(2)-o(1)+v(2,sp) ... 
      -v(1,sp)-w)/10:o(2)+v(2,sp)]);
    set(ss(3,sp),'XData',[o(2)+v(2,sp)+w:(o(3)-o(2)+v(3,sp) ...
      -v(2,sp)-w)/10:o(3)+v(3,sp)]);
  end
  drawnow;                               % draw updated positions
end
