function [x t u] = wave(tmax, level, lambda, al, x0l, deltal, ...
                        ar, x0r, deltar, trace)
% wave Solves 1+1 wave equation using O(h^2) FDA [PHYS210-fda]
%
% wave solves the 1+1 wave equation
%
%     u(t, x)_tt = u_xx
%
% on the domain
%
%     x <= 0 <= 1
%     0 <= t <= tmax
%
% subject to the initial conditions
%
%     u(0, x)   = u_0(x)
%     u_t(0, x) = v_0(x)
% 
% and the (Dirichlet / fixed) boundary conditions 
%
%     u(t, 0) = 0
%     u(t, 1) = 0
%
% The "usual" differencing scheme that employs the 3-pt centred
% formula for the second derivative is used.
% 
% The initial conditions are actually specified in terms of initially
% left- and right-moving components, u_l(x) and u_r(x), respectively,
% each of which is a Gaussian.

   % Internal (debug) tracing.
   ltrace = 1;

   if ltrace
      fprintf('wave: tmax=%g\n', tmax);
      fprintf('wave: level=%d\n', level);
      fprintf('wave: lambda=%g\n', lambda);
      fprintf('wave: al=%g\n', al);
      fprintf('wave: x0l=%g\n', x0l);
      fprintf('wave: deltal=%g\n', deltal);
      fprintf('wave: ar=%g\n', ar);
      fprintf('wave: x0r=%g\n', x0r);
      fprintf('wave: deltar=%g\n\n', deltar);
   end

   err = 0;

   % Define 0-length defaults for output vectors/matrices
   x = zeros(0,0);
   t = zeros(0,0);
   u = zeros(0,0);

   % Check arguments for validity
   if tmax <= 0 
      errmsg('tmax must be > 0'); err = 1;
   end
   if level ~= fix(level)
      fprintf('level is not an integer\n'); err = 1;
   else
      if level < 3 
         fprintf('level must be >= 3\n'); err = 1;
      end
   end
   if lambda <= 0
      fprintf('lambda must be >= 0\n'); err = 1;
   end
   if (x0l < 0.0) | (x0l > 1.0) 
      fprintf('x0l must be in range 0.0 - 1.0\n'); err = 1;
   end
   if (x0r < 0.0) | (x0r > 1.0) 
      fprintf('x0r must be in range 0.0 - 1.0\n'); err = 1;
   end
   % Bail out if there were any errors in the arguments
   if err
      return
   end

   % Define x and t vectors and related mesh parameters
   nx = 2^level + 1;
   x = linspace(0.0, 1.0, nx);
   delx = x(2) - x(1);

   delt1 = lambda * delx;
   nt = fix(tmax / delt1 + 1);
   t = linspace(0.0, tmax, nt);
   delt = t(2) - t(1);

   % Initialize solution array
   u = zeros(nx, nt);

   if ltrace
      fprintf('wave: nx=%d\n', nx);
      fprintf('wave: delx=%g\n', delx);
      fprintf('wave: nt=%d\n', nt);
      fprintf('wave: delt1=%g\n', delt1);
      fprintf('wave: delt=%g\n', delt);
      fprintf('wave: u [%d %d]\n', size(u));
      fprintf('wave: x [%d]\n', length(x));
      fprintf('wave: t [%d]\n', length(t));
   end

   % Initialize first and second time steps using instructor-supplied
   % functions
   u(:, 1) = waveu1(x, delt, al, x0l, deltal, ar, x0r, deltar)';
   u(:, 2) = waveu2(x, delt, al, x0l, deltal, ar, x0r, deltar)';

   % Advance the discrete solution the required number of time steps
   % using the difference equations
   lambdasq = lambda^2;
   if trace
      fprintf('wave: Number of time steps to take: %d\n\n', nt);
      fprintf('wave: Step ');
   end
   for n = 2 : nt-1
      if trace && rem(n, trace) == 0 
         fprintf('  %d', n);
      end
      u(2:nx-1, n+1) = 2 * u(2:nx-1, n) - u(2:nx-1, n-1) + ...
         lambdasq * (u(1:nx-2, n) - 2 * u(2:nx-1, n) + u(3:nx, n));
   end
   if trace
      fprintf('  Done.\n');
   end
end
