% The wave equation in 1D. The PDE is represented as 2N coupled ODEs
% on an N-element spatial grid with periodic boundaries.
function dqdt = waveq(q, t, Delta, c)
% q     = state vector. Should be a row vector, the first half being
%         the wavefunction u, and the second half the time derivative of u.
% t     = time. Not used. Provided for compatibility with ODE solvers.
% Delta = grid spacing.
% c     = wave speed.

N = numel(q)/2; % Number of grid points

% We implement periodic boundary conditions on u by wrapping around both 
% ends, so u has N+2 elements.
u = [ q(N), q(1:N), q(1) ]; 
udot = q( (N+1):(2*N) );

% Here is the PDE:
udotdot = (c/Delta)^2 * ( u(1:N) - 2*u(2:(N+1)) + u(3:(N+2)) );
   % Note how those periodic BCs on u work --- now udotdot has N elements.

dqdt = [udot, udotdot];

