function Prob21p17

% First step is to change 2 second order ODEs into 4 first order ODEs
% Let
% dx1/dt = x3
% dx2/dt = x4
% dx3/dt = d^2x1/dt^2 = first dynamics equation
% dx4/dt = d^2x2/dt^2 = second dynamics equation
% The state vector x is then
% [x1 x2 x3 x4] where
% x1 and x2 are the initial positions and x3 and x4 are the initial
% velocities

k1=5;k2=5;m1=2;m2=2;w1=5;w2=5;L1=2;L2=2;
[t,x]=ode45(@TwoMasses,[0 20],[2 15 0 0],[],k1,k2,m1,m2,w1,w2,L1,L2);
subplot(2,1,1);plot(t,x(:,1),t,x(:,2),'--')
grid;title('displacements');legend('x1','x2')
subplot(2,1,2);plot(t,x(:,3),t,x(:,4),'--')
grid;title('velocities');legend('v1','v2')
pause
subplot(1,1,1),plot(x(:,2),x(:,1))
xlabel('x2');ylabel('x1')

function dx = TwoMasses(t,x,k1,k2,m1,m2,w1,w2,L1,L2)
dx = [x(3);
      x(4);
      -k1/m1*(x(1)-L1)+k2/m1*(x(2)-x(1)-w1-L2);
      -k2/m2*(x(2)-x(1)-w1-L2)];