%What happens when two masses collide?
%Parameters
v1b = 2;   v2b = -2;  % velocities of masses before collision
m1  = 1;   m2  = 1;  % masses
e   = -1;

%Governing Equations (unknowns on left, knowns on right)
%Mom. Cons.:  m1*v1a + m2*v2a  = m1*v1b + m2*v2b
%Restitution;   -v1a +    v2a  = e*(v1b - v2b)

A = [ m1   m2;   %A matrix of coefficients of the ...
      -1    1];  %... unknowns v1a and v2a
  
rhs = [m1*v1b + m2*v2b;  %All of the before stuff that is known
       e*(v1b - v2b)];

% =>         A*va      = rhs     % Matrix form 

%Solve simultaneous equations
va = A\rhs;  %Calculation: Matlab's great BACKSLASH command.

v1a = va(1); v2a=va(2);

Eb = (m1*v1b^2 + m2*v2b^2)/2;
Ea = (m1*v1a^2 + m2*v2a^2)/2;
vcmb = (m1*v1b + m2*v2b)/(m1 + m2);
vcma = (m1*v1a + m2*v2a)/(m1 + m2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(['Coefficient of restitution e =  '            num2str(e)])
disp(['v1b =    ', num2str(v1b) '   v1a =    '  num2str(v1a)  ])
disp(['v2b =    ', num2str(v2b) '   v2a =    '  num2str(v2a)  ])
disp(['Eb  =    ', num2str(Eb)  '   Ea  =    '  num2str(Ea)   ])
disp(['vcmb=    ', num2str(vcmb)  '  vcma  =    '  num2str(vcma)   ])