function exercise3 ( )

%*****************************************************************************80
%
%% exercise3 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, 'exercise3:\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 = nu * uxxp + up .* ( 1.0 -  up );
      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
function lp = laplacian_circle ( up, dx )

%*****************************************************************************80
%
%% laplacian_circle approximates the laplacian on a periodic 1D domain.
%
%  Discussion:
%
%    If the domain is represented by nx equally spaced points, it is
%    assumed that only the data for the first nxp=nx-1 points is supplied.
%    In other words, the last entry of up() is at the node just BEFORE
%    the node where the periodic condition takes effect.
%
%    This is the standard laplacian estimate on a uniform grid, 
%    except that, by periodicity, and suppressing the last node, we have
%    the east neighbor of node 1 is the last node, and
%    the west neighbor of the last node is node 1.
%
%    This function offers four equivalent ways of computing the Laplacian.
%
%  Modified:
%
%    27 June 2020
%
%  Input:
%
%    real up(nxp): the values of a function on an equally spaced grid
%    of points in a periodic domain.
%
%    real dx: the spacing between any two consecutive nodes.
%
%  Output:
%
%    real lp(nxp): the estimate for the laplacian.
%
  nxp = length ( up );
  lp = zeros ( size ( up ) );
%
%  Simplest formula.
%
  if ( false )

    ip = 1;
    lp(ip) = ( up(ip+1) - 2.0 * up(ip) + up(nxp) ) / dx^2;

    for ip = 2 : nxp - 1
      lp(ip) = ( up(ip+1) - 2.0 * up(ip) + up(ip-1) ) / dx^2;
    end

    ip = nxp; 
    lp(ip) = ( up(1) - 2.0 * up(ip) + up(ip-1) ) / dx^2;
%
%  Use the mod() function to take care of indexing.
%
  elseif ( false )

    for ip = 1 : nxp
      east = 1 + mod ( ip,     nxp );
      west = 1 + mod ( ip - 2, nxp );
      lp(ip) = ( up(east) - 2.0 * up(ip) + up(west) ) / dx^2;
    end
%
%  Use "circshift" to force indices to wrap around.
%
  elseif ( true )
    lp = (         circshift ( up,  1 ) ... 
           - 2.0 * circshift ( up,  0 ) ...
                 + circshift ( up, -1 ) ) / dx^2;
%
%  Just make a list of the indexes.
%
  else

    im1 = [ nxp, 1:nxp-2, nxp-1 ];
    i =   [ 1,   2:nxp-1, nxp   ];
    ip1 = [ 2,   3:nxp,   1     ];
    lp = ( up(ip1) - 2.0 * up(i) + up(im1) ) / dx^2;

  end

  return
end

