% Quantum particle in a box.
% The time dependent Schrödinger equation is finite differenced in space to
% become an array of N coupled ODEs. Any ODE solver can then be used to
% integrate in time.
function dqdt = particle_box(q, t)
% q is the wave function (state vector).
% t is dimensionless time.
% The domain is assumed to run from x=0 to x=L, where L=1 is hard coded.
%
% The PDE implemented here is: 
%   dq/dt = 0  at the ends (Dirichlet BC), and
%   dq/dt = i d^2 q / dx  in the interior.
% The phase speed is c = k, where k is the wavenumber in radians per 
% unit length. [The units appear mismatched because of the dimensionless
% units of the PDE here.]
% The maximum wavenumber is the Nyquist frequency,
%   k_max = pi * N / L = N pi = c_max.
% So, the Courant condition tells us that the longest usable timestep is
%   dt_max = dx/c_max = (L/(N-1))/(N*pi) ~ 1 / (N^2 * pi).
%
% By contrast, the slowest mode would have period L and wavenumber
%   k_min = c_min = 2*pi/L = 2*pi.
% The crossing time for the slowest mode is then
%   t_max = L/c_min = 1/(2*pi).
%
% So, a simulation allowing a full period of the slowest mode would have to 
% take a large number of steps:
%   N_steps = t_max / dt_max = N^2 / 2.

N = numel(q); % resolution of domain, which is up to the person setting up the problem.
L = 1.0       % length of domain
dx = L/(N-1); % uniform sampling pitch

% Ensure that q satisfies the homogeneous Dirichlet boundary conditions.
q(1) = 0;
q(N) = 0;

% Schrodinger equation for zero potential, with above BCs.
dqdt = i * [0, q(1:N-2) - 2*q(2:N-1) + q(3:N), 0]/(dx.^2);
