% Telegrapher's equations for sound in 1D. PDE is represented as 2N coupled ODEs
% on an N-element spatial grid with periodic boundaries.
function dqdt = telegrapher(q, t, Delta, c)
% q     = state vector. Should be a row vector, the first half being
%         the dimensionless velocity v/c, and the second half the 
%         dimensionless density (rho-<rho>)/<rho>.
% t     = time. Not used. Provided for compatibility with ODE solvers.
% Delta = grid spacing.
% c     = wave speed.
%
% To set up initial conditions: velocity = density gives a wave traveling
% in the +x direction, while velocity = -density gives a wave traveling in the
% -x direction.

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

% We implement periodic boundary conditions by wrapping around both 
% ends, so each of our two variables has N+2 elements.
velocity = [ q(N), q(1:N), q(1) ]; 
density = [ q(2*N), q((N+1):(2*N)), q(N+1) ];

% PDEs:
dt_vel = - c * (  density(3:(N+2)) -  density(1:N) ) / (2*Delta);
dt_den = - c * ( velocity(3:(N+2)) - velocity(1:N) ) / (2*Delta);
   % Note how those periodic BCs on u work --- time derivs have N elements.

dqdt = [dt_vel, dt_den];

