% Wave equation in a 1D box.
% The wave 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 = waveq_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 state of the system comprises z(x,t) and v(x,t).
% The equations of motion and BCs are: 
%   z = v = 0           at the ends (Dirichlet BC).
%   dz/dt = v           everywhere.
%   dv/dt = d^2 q / dx  in the interior.
%   dv/dt = 0           at the ends.
% The phase speed is c = 1 (z, v, x, t are nondimensional)
% The Courant condition tells us that the longest usable timestep is
%   dt_max = dx/c = (L/(N-1)) ~ 1/N.
% The crossing time, which is necessarily also the period of the fundamental mode, is
%   t_max = L/c = 1.
% So, a simulation allowing a full period of the lowest frequency mode would have to 
% take the following number of steps:
%   N_steps = t_max / dt_max = N.

N = numel(q)/2; % 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

%Extract wave function and its rate of change from q.
z = q(1:N);
v = q(N+1:2*N);

% Ensure that z and v satisfy the homogeneous Dirichlet boundary conditions.
z(1) = 0;
z(N) = 0;
v(1) = 0;
v(N) = 0;

% Wave equation, with above BCs.
dzdt = v;
dvdt = [0, z(1:N-2) - 2*z(2:N-1) + z(3:N), 0]/(dx^2);
dqdt = [dzdt, dvdt];
