function advection_pde_ftcs ( )

%*****************************************************************************80
%
%% advection_pde_ftcs() solves the advection PDE using the FTCS method.
%
%  Discussion:
%
%    The FTCS method (Forward Time, Center Space) is unstable for the 
%    advection problem.
%
%    Given a smooth initial condition, successive FTCS approximations will
%    exhibit erroneous oscillations of increasing magnitude.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 May 2021
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Willem Hundsdorfer, Jan Verwer,
%    Numerical solution of time-dependent advection-diffusion-reaction equations,
%    Springer, 2003
%    ISBN: 978-3-662-09017-6
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'advection_pde_ftcs():\n' );
  fprintf ( 1, '  Solve the constant-velocity advection PDE in 1D,\n' );
  fprintf ( 1, '    du/dt + c du/dx = 0\n' );
  fprintf ( 1, '  over the interval:\n' );
  fprintf ( 1, '    0.0 <= x <= 1.0\n' );
  fprintf ( 1, '  with periodic boundary conditions:\n' );
  fprintf ( 1, '    u(0) = u(1)\n' );
  fprintf ( 1, '  and advection velocity\n' );
  fprintf ( 1, '    c = constant\n' );
  fprintf ( 1, '  and initial condition\n' );
  fprintf ( 1, '    u(0,x) = (10x-6)^2 (8-10x)^2 for 0.6 <= x <= 0.8\n' );
  fprintf ( 1, '           = 0 elsewhere.\n' );
  fprintf ( 1, '  and NX equally spaced nodes in X,\n' );
  fprintf ( 1, '  and NT equally spaced points in T,\n' );
  fprintf ( 1, '  using the FTCS method:\n' );
  fprintf ( 1, '    FT: Forward Time  : du/dt = (u(t+dt,x)-u(t,x))/dt\n' );
  fprintf ( 1, '    CS: Centered Space: du/dx = (u(t,x+dx)-u(t,x-dx))/2/dx\n' );
%
%  Get the parameters.
%
  [ c, t0, tstop, xmin, xmax ] = advection_parameters ( );

  nx = 101;
  dx = ( xmax - xmin ) / ( nx - 1 );
  x = linspace ( xmin, xmax, nx );

  nt = 8001;
  dt = ( tstop - t0 ) / ( nt - 1 );
  t = linspace ( t0, tstop, nt );

  h = zeros ( nt, 1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes NX = %d\n', nx );
  fprintf ( 1, '  Number of time steps NT = %d\n', nt );
%
%  Setting up im1, i, ip1 allows us to define the centered spatial difference
%  with the periodic boundary condition.
%
  Im1 = [ nx, 1:nx-2, nx-1 ];
  I   = [ 1,  2:nx-1, nx   ];
  Ip1 = [ 2,  3:nx,   1    ];
%
%  Compute the solution at each time.
%
  for j = 1 : nt
%
%  Apply the initial condition to get solution at initial time.
%
    if ( j == 1 )
      u = advection_initial ( x );
      umax = max ( u ) * ones ( 1, nx );
      umin = min ( u ) * ones ( 1, nx );
%
%  At later times, rearrange
%
%    u(x(i),t(j)) - u(x(i),t(j-1))       u(x(i+1),t(j-1)) - u(x(i-1),t(j-1))
%    ----------------------------- + c * ----------------------------------- = 0
%                dt                                    2 dx
%
%  to get
%
%    u(x(i),t(j)) = u(x(i),t(j-1)) - c * dt / ( 2 dx ) * ( u(x(i+1),t(j-1)) - u(x(i-1),t(j-1)) )
%
    else
      u(I) = u(I) - c * dt / dx / 2.0 * ( u(Ip1) - u(Im1) );
    end

    h(j) = advection_conserved ( x, u );
%
%  Display the solution.
%
    plot ( x, umin, 'r--', ...
           x, umax, 'r--', ...
           x, u, 'b-', ...
           'linewidth', 3 );

    grid ( 'on' );
    xlabel ( '<--- x --->' );
    ylabel ( '<--- u(x,t) --->' );
    axis ( [ 0.0, 1.0, -2.0, +2.0 ] )
    title_string = sprintf ( 'Advection PDE FTCS: Time = %g\n', t(j) );
    title ( title_string );

    if ( mod ( j, 800 ) == 1 )
      filename = sprintf ( 'advection_pde_ftcs_%d.png', j );
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    end

  end
%
%  Plot the conservation quantity.
%
  figure ( );
  clf ( );
  hold ( 'on' );
    plot ( t, h, 'r-',  'linewidth', 3 );
    plot ( [ t0, tstop ], [ h(1), h(1) ], 'k:',  'linewidth', 3 );
    plot ( [ t0, tstop ], [ 0.0, 0.0 ], 'k--',  'linewidth', 3 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<--- t --->' );
  ylabel ( '<--- h(t) --->' );
  legend ( 'H computed', 'H initial', 'X axis', 'Location', 'South' );
  title ( 'advection pde FTCS: conservation' );
  filename = 'advection_pde_ftcs_conservation.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end

