Visualizing the (expected) gradual change from quantum to classical behavior I'm pretty much a layman (no formal Physics education beyond high school) but find Quantum Mechanics endlessly fascinating.
Ok, so when the double slit experiment is performed using photons/electrons/etc we get a proper interference pattern.
If it's performed using much larger objects, say footballs (which also have dual nature), we'd get (in 99.9999999...% cases) no interference pattern — just two bands corresponding to the two slits.
Now, one would expect that as the size of the particles is gradually increased, there'd be a gradual change from the pattern created using photons/electrons (with a maxima band right behind the middle of the slits, and smaller alternating bands, etc) to  the pattern for large/classical objects (with just two bands corresponding to the two slits). 
Is there a way to see this change happening in an animation? If I wanted to create one programmatically, are there good frameworks/libraries to easily do this in some language?
Edit:
Bonus points if you can make the animation. Or, at the least, provide a nice verbal description of how this transition would happen.
 A: I love making small pieces of simulation. I would recommend octave/Matlab for such a simple simulations. To make that point, here is a small piece of code I wrote in 10 minutes with octave/Matlab. It simulates a double slit experiment by solving a 2D Schrödinger equation in a box with a "double-slit source term" using finite difference and Crank-Nicholson propagator.
In principle, it plots the solution to this partial differential equation (except that I forgot do multiply the source term with dt in the code, but that doesn't matter here) 
$$ i\partial_t \Psi(t) = (-\frac{1}{2}\nabla^2) \Psi(t) + e^{iwt}(\delta(r-r_{slit1}) + \delta(r-r_{slit2})), \Psi(0)=0.$$
It uses atomic units, where electron mass is assumed to be 1. If you would like to see how going classical affects, you would have to divide with mass (in units of electron mass) in definition of kinetic energy operator L.
To avoid unwanted reflections, the source terms are in the middle of the box. Also, the detector is far box box walls to avoid reflections. The code could be improved by adding absorbing boundary conditions.
% Define a 2D grid
g = 150; 
p=linspace(-10,10,g);
h=p(2)-p(1)

% Define laplace operator in 1D
e=ones(g,1);
L=spdiags([e -2*e e], -1:1, g, g)/h^2;

% Create 2D kinetic energy operator
I = speye(g);
I2 = kron(I,I);
T = -0.5*(kron(L,I) + kron(I,L));

dt = 0.02;

% Initial wave function 
psi0=zeros(g,g);

w = 8; % Frequency of the wave
t=0;

A1=1;
A2=1;
while 1

% Real part of the wave function
% TODO: fixed bounds
subplot(2,2,1);
imagesc(real(psi0));
line([40 40],[1 g]);
title('Real wf');

% Imaginary part of the wave function
% TODO: fixed bounds
subplot(2,2,2);
imagesc(imag(psi0));
line([40 40],[1 g]);
title('Imaginary wf');

% Simulated measurement device
subplot(2,2,3);
plot(abs(psi0(:, 40)).^2);
title('Intensity at measurement device');

% Add source terms for the double slit
psi0(g/2-8,g/2)=psi0(g/2-8,g/2) + A1*exp(1i*w*t);
psi0(g/2+8,g/2)=psi0(g/2+8,g/2) + A2*exp(1i*w*t);

% Change from matrix to vector to apply propagation
psi0=psi0(:);
% Propagate wave function using Crank-Nicholson propagation
% Equal to psi0=expm(-1i*T*dt)*psi0, but faster
psi0=(I2-1i*T*dt/2)\((I2+1i*T*dt/2)*psi0);
% Change back to 2D matrix
psi0=reshape(psi0,g,g);

t=t+dt;
pause(0.001);
end


