function diffusion_pde_ftcs ( )

%*****************************************************************************80
%
%% diffusion_pde_ftcs() solves the diffusion PDE using the FTCS method.
%
%  Discussion:
%
%    The FTCS method (Forward Time, Center Space) is used.
%
%  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 diffusion-diffusion-reaction equations,
%    Springer, 2003
%    ISBN: 978-3-662-09017-6
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'diffusion_pde_ftcs():\n' );
  fprintf ( 1, '  Solve the diffusion PDE in 1D,\n' );
  fprintf ( 1, '    du/dt - mu d2u/dx2 = 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 diffusion coefficient\n' );
  fprintf ( 1, '    mu = 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: d2u/dx2 = (u(t,x+dx)-2u(t,x)+u(t,x-dx))/dx^2\n' );

  [ mu, t0, tstop, xmin, xmax ] = diffusion_parameters ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Parameter values:\n' );
  fprintf ( 1, '    mu    = %g\n', mu );
  fprintf ( 1, '    t0    = %g\n', t0 );
  fprintf ( 1, '    tstop = %g\n', tstop );
  fprintf ( 1, '    xmin  = %g\n', xmin );
  fprintf ( 1, '    xmax  = %g\n', xmax );

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

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

  h = zeros ( nt, 1 );

  cfl = mu * dt / dx^2;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes NX = %d\n', nx );
  fprintf ( 1, '  Number of time steps NT = %d\n', nt );
  fprintf ( 1, '  CFL coefficient ( must be < 0.5 ) = %g\n', cfl );
%
%  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    ];
%
%  Initialize the movie.
%
  colormap ( 'jet' );
  moviename = 'diffusion_pde_ftcs.avi';
  Video = VideoWriter ( moviename, 'Motion JPEG AVI' );
  Video.FrameRate = 10;
  open ( Video );
%
%  Compute the solution at each time.
%
  for j = 1 : nt
%
%  Apply the initial condition to get solution at initial time.
%
    if ( j == 1 )
      u = diffusion_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)) - 2u(x(i),t(j-1)) + u(x(i-1),t(j-1))
%    ----------------------------- - mu * ----------------------------------------------------- = 0
%                dt                                               dx^2
%
%  to get
%
%    u(x(i),t(j)) = u(x(i),t(j-1)) ...
%      + mu * dt / ( 2 dx^2 ) * ( u(x(i+1),t(j-1)) - 2 u(x(i),t(j-1) + u(x(i-1),t(j-1)) )
%
    else
      u(i) = u(i) + mu * dt / dx^2 * ( u(ip1) - 2 * u(i) + u(im1) );
    end

    h(j) = diffusion_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, -0.05, 1.05 ] )
    title_string = sprintf ( 'diffusion PDE FTCS, Time = %8.4f\n', t(j) );
    title ( title_string );

    if ( j == 1 )
      filename = 'diffusion_pde_ftcs_initial.png';
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    elseif ( j == nt )
      filename = 'diffusion_pde_ftcs_final.png';
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    end

    frame = getframe ( 1 );
    writeVideo ( Video, frame );

    pause ( 0.125 );

  end
%
%  Close the movie.
%
  close ( Video );
  fprintf ( 1, '  Movie saved as "%s"\n', moviename );
%
%  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 ], 'k-', 'linewidth', 3 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<--- t --->' );
  ylabel ( '<--- h(t) --->' );
  legend ( 'H computed', 'H initial', 'X axis', 'Location', 'South' );
  title ( 'diffusion PDE ftcs: conservation' );
  filename = 'diffusion_pde_ftcs_conservation.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end

