% this program plots the displacement of each mass. positive displacement
% means that the objects is moved to the right, negative means to the left.


nsprings=20;  % number of masses (numnber of springs is one less
             % than the number of masses)
x=1:nsprings+1;   % initial position of all of the masses.

vx=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
a=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
% initial velocities are all zero and so are accelerations

k=[.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1,.1];
% spring constant of each spring
% spring "0" is the spring connecting mass(1) [the wall] and mass(2),
%  and so on.
% masses are all assumed to be 1.

dt=0.01;   % is this a good time step?
x(2)=2.5;  % <=== here is the initial condition!
%initial condition is that one mass is stretched to a new position.

displacement=x-(1:21)    % calculate displacement from initial positions
h=plot(displacement);
xlabel('mass number')
ylabel('displacement')
axis([1 21 -1 1])
set(h,'EraseMode','none','LineWidth',3)


while 1  		% run the simulation until you hit cntrl-C
   % first, calculate displacement of each spring on the left
   %  of each mass
   stretch=(x(2:nsprings)-x(1:nsprings-1))-1;
   % and the acceleration from Hooke's law
   a(2:nsprings)=-k.*stretch;    % note the matrix multiply!
   % but the wall is fixed, so there can't be any acceleration at the wall.
   a(nsprings)=0;
  
   % now look at the spring to the right of each mass and add that acceleration
   stretch=(x(1:nsprings-1)-x(2:nsprings))+1;
   a(1:nsprings-1)=a(1:nsprings-1)-k.*stretch;
   a(1)=0;	% the left wall doesn't move either

   oldx=x;
   % calculate new positions and velocities
   x=.5*a.*dt^2+vx.*dt+x;		% a little simple physics here!
   vx=vx+a.*dt;


   displacement=x-(1:21);
   set(h,'YData',displacement,'Color',[0,0,1])  % overplot a blue line
   set(h,'Color',[1,0,0])               % overplot a red line
   drawnow                              % plot it!
end