function [U,M,flag]=wavesolve(c,tmax,nt,nx,u0,ut0,myf)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% This program solves the problem 
%   u_{tt} - c^2 u_{xx} = f(x,t)
% where 0 < x < 1 and t 0 <= t <= tmax.  The boundary conditions are
% u(0,0) = u(1,0) = 0.
% 
% On input:
%    c > 0 is the speed of propagation of the wave 
%      tmax (integer) is the upper bound on the values of t of interest
%      nt is the (integer) number of timesteps per unit time
%      nx is the (integer) number of steps per unit x 
%    u0 is a vector of function values for u(x_k,0),
%    ut0 is a vector of function values for u_t(x_k,0).
%     with $x_k = k deltax$ and k = 1,...,nx-1.
%     deltax = 1/nx  is the stepsize in x.
%     deltat = 1/nt  is the stepsize in t.
%    f is a pointer to a function that evaluates f(x,t).
% On output:
%     U is an array, with element (k,j) giving the computed value
%           of u(k deltax, (j-1)deltat)
%             k = 1,...,nx-1    j = 0,... tmax/deltat.
%     flag = 0 if no error occurred. 
%          = 1 if  c, tmax, nt, or nx <= 0
%              or tmax, nt, or nx not an integer
%          = 2 if stability condition c nx / nt <= 1 violated.
%     M is a Matlab movie displaying the solution.
% 
% The program uses 2nd order finite differences in space and time.
% wavesolve.m     Dianne P. O'Leary  05/2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Check validity of input and set the error flag.

flag = 0;
if (c <= 0)
  disp(sprintf('wavesolve c = %e should be positive.',c))
  flag = 1;
end
if (tmax <= 0) | (floor(tmax) ~= tmax)
  disp(sprintf('wavesolve tmax = %e should be a positive integer.',tmax))
  flag = 1;
end
if (nt <= 0) | (floor(nt) ~= nt)
  disp(sprintf('wavesolve nt = %e should be a positive integer.',nt))
  flag = 1;
end
if (nx <= 0) | (floor(nx) ~= nx)
  disp(sprintf('wavesolve nx = %e should be a positive integer.',nx))
  flag = 1;
end

% Initialize some constants and check the stability condition.

deltat = 1/nt;
deltat2 = deltat^2;
deltax = 1/nx;
deltax2 = deltax^2;
if (c*deltat/deltax  >=  1) 
  disp(sprintf('wavesolve stability condition violated.'))
  disp(sprintf('  Increase the ratio of nx to nt.      '))
  flag = 2;
end
c2 = c*c;
xext = linspace(0,1,nx+1)';

nstep = tmax*nt;

% Initialize U to keep Matlab from repeatedly allocating space.

U = zeros(length(u0),nstep+1);

% The iteration is started using a 2nd order finite difference
% approximation to u_t(k deltax,0):
% (u(k deltax, deltat) - u(k deltax, -deltat))/(2 deltat) 
%                    approx= u_t(k deltax, 0).
% This gives us an expression for u(k deltax, -deltat) that
% we use to compute u(k deltax, (j+1) deltat) in the iteration
%
% [u(k deltax, (j-1)deltat) - 2 u(k deltax, j deltat) 
%              + u(k deltax, (j+1)deltat)] / deltat^2 
%  - c^2 [u((k-1)deltax, j deltat) -2 u(k deltax, j deltat) 
%              + u((k+1)deltax, j deltat) ] / deltax^2 
%   = f(k delta x, j delta t)

U(:,1) = u0;
ustep = 2*u0 + 2*deltat*ut0;
fnow = f(xext,0);
fval = fnow(2:end-1);
uext = [0;U(:,1);0];
udiff = c2*(-uext(3:end) + 2 *uext(2:end-1) - uext(1:end-2))/deltax2;
U(:,2) = (ustep + deltat2*(fval - udiff))/2;

% Continue the iteration, now that we have values at two previous
% timesteps.

for j=3:nstep+1,

  ustep = 2*U(:,j-1) - U(:,j-2);
  fnow = f(xext,(j-2)*deltat); % evaluate at previous timestep
  fval = fnow(2:end-1);
  uext = [0;U(:,j-1);0];
  udiff = c2*(-uext(3:end) + 2 *uext(2:end-1) - uext(1:end-2))/deltax2;
  U(:,j) = ustep + deltat2*(fval - udiff);

end

% Find the max and min values of U to scale the plots.

umax = max(max(U));
umin = min(min(U));

% Display the plots and make the movie.

for j=1:nstep+1
   plot(xext,[0;U(:,j);0])
   axis([0 1 umin umax])
   title('dummy')
   M(j) = getframe;
end

% Three test problems were to be run.
%
% Struck string:  The true solution is a square wave that travels to the
% right, reflects to its negative, travels to the left, etc.
% (Picture: p. 539, Kammler)
%
% Plucked string: The v-shaped wave created by the pluck travels to the right,
% diminishing in amplitude.  Meanwhile, a growing (but negative) v-shaped wave
% follows it.  The process then repeats in reverse, with the negative
% wave traveling left and diminishing while a growing positive one follows,
% and so on.
% (similar Picture: p. 539, Kammler)
%
% Forced vibration: We excited an eigenmode, so we see the initial condition
% persist but grow in amplitude.
% There is a node at x=.5