function skeleton3 ( )

%*****************************************************************************80
%
%% skeleton3 solves a version of Fisher's equation.
%
%  Discussion:
%
%    The equation is dudt = nu * duxx + u*(1-u).  
%
%    The initial values should be between 0 and 1.  
%
%    Any nonzero values quickly tend to u=1.
%
%    Periodic boundary conditions are used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'skeleton3:\n' );
  fprintf ( 1, '  Solve a version of Fisher''s reaction-diffusion equation.\n' );
%
%  Diffusion coefficient (must be nonnegative)
%  Compare the solution behavior for nu = 0, nu = 0.001, and nu = 1.
%
  nu = 1.0;
  nu = 0.001;
%
%  Spatial discretization.
%
  nx = 76;
  xmin = 0.0;
  xmax = 75.0;
  x = linspace ( xmin, xmax, nx );
  dx = ( xmax - xmin ) / ( nx - 1 );
%
%  Time discretization.
%
  nt = 61;
  tmin = 0.0;
  tmax = 15.0;
  t = linspace ( tmin, tmax, nt );
  dt = ( tmax - tmin ) / ( nt - 1 );
%
%  Set up the periodic domain by dropping the last point.
%
  nxp = nx - 1;
  xp = x(1:nxp);
%
%  The time loop.
%
  for k = 1 : nt
%
%  Initial condition involves two step functions.
%
    if ( k == 1 )

      up = 1.00 * ( xp <= 15 ) ...
         + 0.25 * ( 35.0 <= xp & xp <= 45.0 );
%
%  Take a forward Euler time step:
%    dudt = nu duxx + u*(1-u).
%
    else
      uxxp = laplacian_circle ( up, dx );
      dudt = ?
      up = up + dt * dudt;
    end
%
%  For plotting, use the full domain (x,u).
%
    u(1:nxp) = up;
    u(nx) = up(1);

    plot ( x, u, 'b-', 'linewidth', 3 );
    grid ( 'on' );
    xlabel ( '<--- x --->' );
    ylabel ( '<--- u(x,t) -->' );
    label = sprintf ( 'Fisher: k = %3d, t = %4.1f', k, t(k) );
    title ( label );
    axis ( [ xmin, xmax, 0.0, 1.0 ] );
    pause ( 0.1 );

  end

  filename = 'exercise3.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end

