I recently developed a simulation in Matlab of two molecules colliding under the van der waals potential so students could study how elastic collisions worked. The code for the calculation is at the bottom of the page.

Students were given different parameters to try and test if momentum and kinetic energy were conserved. Here is the general feedback from the lab:

• Things can collide without touching.
• Many students didn’t do any calculations. They just said it was conserved and moved on. 🙁
• Many students forgot that momentum is a vector and you need to test the x- and y-directions independently.
• On the flip side, many students also forgot that Kinetic Energy is a scalar. Therefore it doesn’t have components. Find the magnitude of the velocity by calculating. v² = v_x² + v_y² .

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Models 2-D elastic van-der Walls Collisions between Molecules    %%%%
%%% By Matthew Wright                                                %%%%
%%%  U(r) = c6 / r^6 + c12 / r^12                                    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all
clear all

%% Parameters for the problem
x1    = -10;  %x-Position 1
y1    = 0;    %y-Position 1
x2    = 0;    %x-Position 2
y2    = 8;    %y-Position 2
v1x   = 10;    %x-velocity 1
v1y   = 0;    %y-velocity 1
v2x   = 0;    %x-velocity 2
v2y   = 0;    %y-velocity 2
m1    = 1;    %mass 1 
m2    = 1;    %mass 2
dt    = 0.01; %time step
c6    = -1e2;    %Van-Der Walls Parameter
c12   = +1e8;    %Van-Der Walls Parameter

%% Build the picture
plot ([-10, 10], [-10, 10], 'w')
rec1 = rectangle ('Position', [x1,y1,1,1], 'FaceColor', 'r', 'Curvature', [1,1]);
rec2 = rectangle ('Position', [x2,y2,1,1], 'FaceColor', 'b', 'Curvature', [1,1]);
axis ([-10 10 -10 10]);

time = 0;
xx1  = [];
yy1  = [];
xx2  = [];
yy2  = [];
while time < 10
    
    %% Get the force
    rr = sqrt((x2-x1)^2+(y2-y1)^2);
    FF = -6*c6/rr^7 - 12*c12/rr^13;
    Fx = FF * (x2-x1) / rr; %Break into components
    Fy = FF * (y2-y1) / rr; %Break into components
    
    %% Get the acceleration
    a1x = Fx / m1;
    a1y = Fy / m1;
    a2x = -1*Fx / m2; %Netwon's 3rd law
    a2y = -1*Fy / m2; %Netwon's 3rd law  
    
    %% Update velocities
    v1x = v1x + a1x * dt;
    v1y = v1y + a1y * dt;
    v2x = v2x + a2x * dt;
    v2y = v2y + a2y * dt;     
    
    %% Update Positions
    x1 = x1 + v1x * dt;
    x2 = x2 + v2x * dt;
    y1 = y1 + v1y * dt;
    y2 = y2 + v2y * dt;    

%     %% Save Positions
%     xx1 = [xx1, x1];
%     xx2 = [xx2, x2];
%     yy1 = [yy1, y1];
%     yy2 = [yy2, y2];
    
    time = time + dt;
    
    %update the position of the rectangles 
    set(rec1,'Position', [x1,y1,1,1]);
    set(rec2,'Position', [x2,y2,1,1]);
    
%     hold on
%     plot (xx1, yy1, '-r');
%     plot (xx2, yy2, '-b');
%     hold off

    drawnow
end
% Force = @(r) -6*c6/r^7 - 12*c12/r^13; 

v1x
v1y
v2x
v2y