function u = fd1d_wave ( x_num, x1, x2, t_num, t1, t2, c, u_x1, u_x2, ...
  u_t1, ut_t1 )

%*****************************************************************************80
%
%% fd1d_wave() computes a finite difference solution of the 1D wave equation.
%
%  Discussion:
%
%    This program solves the 1D wave equation of the form:
%
%      c^2 Uxx = Utt
%
%    over the spatial interval [X1,X2] and time interval [T1,T2],
%    with initial conditions:
%
%      U(T1,X)  = U_T1(X),
%      Ut(T1,X) = UT_T1(X),
%
%    and boundary conditions of Dirichlet type:
%
%      U(T,X1) = U_X1(T),
%      U(T,X2) = U_X2(T).
%
%    The value C represents the propagation speed of waves.
%
%    The program uses the finite difference method, and marches
%    forward in time, solving for all the values of U at the next
%    time step by using the values known at the previous two time steps.
%
%    Central differences may be used to approximate both the time
%    and space derivatives in the original differential equation.
%
%    Thus, assuming we have available the approximated values of U
%    at the current and previous times, we may write a discretized
%    version of the wave equation as follows:
%
%      Uxx(T,X) = ( U(T,   X+dX) - 2 U(T,X) + U(T,   X-dX) ) / dX^2
%      Utt(T,X) = ( U(T+dt,X   ) - 2 U(T,X) + U(T-dt,X   ) ) / dT^2
%
%    If we multiply the first term by C^2 and solve for the single
%    unknown value U(T+dt,X), we have:
%
%      U(T+dT,X) =        (     C^2 * dT^2 / dX^2 ) * U(T,   X+dX)
%                  +  2 * ( 1 - C^2 * dT^2 / dX^2 ) * U(T,   X   )
%                  +      (     C^2 * dT^2 / dX^2 ) * U(T,   X-dX)
%                  -                                  U(T-dT,X   )
%
%    (Equation to advance from time T to time T+dT, except for FIRST step!)
%
%    However, on the very first step, we only have the values of U
%    for the initial time, but not for the previous time step.
%    In that case, we use the initial condition information for dUdT
%    which can be approximated by a central difference that involves
%    U(T+dT,X) and U(T-dT,X):
%
%      dU/dT(T,X) = ( U(T+dT,X) - U(T-dT,X) ) / ( 2 * dT )
%
%    and so we can estimate U(T-dT,X) as
%
%      U(T-dT,X) = U(T+dT,X) - 2 * dT * dU/dT(T,X)
%
%    If we replace the "missing" value of U(T-dT,X) by the known values
%    on the right hand side, we now have U(T+dT,X) on both sides of the
%    equation, so we have to rearrange to get the formula we use
%    for just the first time step:
%
%      U(T+dT,X) =   1/2 * (     C^2 * dT^2 / dX^2 ) * U(T,   X+dX)
%                  +       ( 1 - C^2 * dT^2 / dX^2 ) * U(T,   X   )
%                  + 1/2 * (     C^2 * dT^2 / dX^2 ) * U(T,   X-dX)
%                  +  dT *                         dU/dT(T,   X   )
%
%    (Equation to advance from time T to time T+dT for FIRST step.)
%
%    It should be clear now that the quantity ALPHA = C * DT / DX will affect
%    the stability of the calculation.  If it is greater than 1, then
%    the middle coefficient 1-C^2 DT^2 / DX^2 is negative, and the
%    sum of the magnitudes of the three coefficients becomes unbounded.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    23 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    George Lindfield, John Penny,
%    Numerical Methods Using MATLAB,
%    Second Edition,
%    Prentice Hall, 1999,
%    ISBN: 0-13-012641-1,
%    LC: QA297.P45.
%
%  Input:
%
%    integer X_NUM, the number of intervals in the X direction.
%    The number of nodes is X_NUM + 1.
%
%    real X1, X2, the first and last X coordinates.
%
%    integer T_NUM, the number of time steps.  The total number
%    of times at which a solution is known is T_NUM + 1.
%
%    real T1, T2, the first and last T coordinates.
%
%    real C, a parameter which gives the speed of waves.
%
%  Output:
%
%    real U(T_NUM+1,X_NUM+1), the solution of the discretized 
%    wave equation for the discrete set of positions and times.
%
  t_delta = ( t2 - t1 ) / t_num;
  x_delta = ( x2 - x1 ) / x_num;
  alpha = c * t_delta / x_delta;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Stability condition ALPHA = C * DT / DX = %f\n', alpha );

  if ( 1.0 < abs ( alpha ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FD1D_WAVE - Warning!\n' );
    fprintf ( 1, '  The stability condition |ALPHA| <= 1 fails.\n' );
    fprintf ( 1, '  Computed results are liable to be inaccurate.\n' );
  end

  u = zeros ( t_num + 1, x_num + 1 );
%
%  The boundary conditions.
%
  u(:,1)       = u_x1 ( t_num );
  u(:,x_num+1) = u_x2 ( t_num );
%
%  The initial data.
%
  u(1,:)       = u_t1 ( x_num );
  ut(1,:)      = ut_t1 ( x_num );
%
%  For the first step in time, we cannot refer to two previous
%  spatial rows of solution data.  Instead, we must use the
%  initial derivative information.
%
  t = 1;
  for x = 2 : x_num
    u(t+1,x) =       alpha^2   * u(t,x+1) / 2 ...
           + ( 1.0 - alpha^2 ) * u(t,x) ...
           +         alpha^2   * u(t,x-1) / 2 ...
           +         t_delta   * ut(t,x);
  end
%
%  All subsequent steps in time rely on two previous rows of solution data.
%
  for t = 2 : t_num
    for x = 2 : x_num

      u(t+1,x) =               alpha^2   * u(t,  x+1) ...
               + 2.0 * ( 1.0 - alpha^2 ) * u(t,  x  ) ...
               +               alpha^2   * u(t,  x-1) ...
               -                           u(t-1,x  );

    end
  end

  return
end
