function skeleton5 ( )

%*****************************************************************************80
%
%% skeleton5 solves the 2D Gray-Scott reaction-diffusion equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'skeleton5\n' );
  fprintf ( 1, '  Solve the 2D Gray-Scott reaction-diffusion equation.\n' );
%
%  Set the grid size.
%
  n = 200;
%
%  Use the same grid size for x and y.
%
  nx = n;
  ny = n;
%
%  Assume the grid spacing is 1.
%
  dx = 1.0;
  dy = 1.0;
%
%  Set the diffusion coefficients for U and V.
%
  du = 0.3;
  dv = 0.15;
%
%  There are several interesting choices for the feed and kill rates.
%
%  Rings.
%
  if ( false )
    f = 0.01;
    k = 0.1;
%
%  Splits.
%
  elseif ( false )
    f = 0.0367;
    k = 0.0649;
%
%  Coral.
%
  elseif ( true )
    f = 0.0545;
    k = 0.0620;
%
%  Bubbles.
%
  elseif ( false )
    f = 0.012;
    k = 0.05;
%
%  Waves.
%
  elseif ( false )
    f = 0.025;
    k = 0.05;
  end
%
%  Set the time step.
%
  dt = 1.0;
%
%  Choose the initial condition.
%
  if ( true )
    initial_condition = 'wavefront';
  elseif ( false )
    initial_condition = 'bar';
  else
    initial_condition = 'square';
  end
%
%  Set the number of time steps.
%
  nt = 15001;
%
%  Set the frame rate for the movie.
%
  frameRate = 25;

  for it = 1 : nt

    if ( it == 1 )

      u = ones ( nx, ny );
      v = zeros ( nx, ny );

      xm = floor ( nx / 2 );
      ym = floor ( ny / 2 );

      if ( initial_condition == 'wavefront' )
        u(xm-30:xm+30,ym-6:ym-4) = 0.75;
        u(xm-30:xm+30,ym-3:ym-1) = 0.50;
        u(xm-30:xm+30,ym  :ym+2) = 0.25;
        u(xm-30:xm+30,ym+3:ym+5) = 0.00;
        v(xm-30:xm+30,ym+3:ym+5) = 1.00;
      elseif ( initial_condition == 'bar' ) 
        u(xm-10:xm+10,ym-2:ym+2) = 0.0;
        v(xm-10:xm+10,ym-2:ym+2) = 1.0;
      else 
        u(xm-5:xm+5,ym-5:ym+5) = 0.0;
        v(xm-5:xm+5,ym-5:ym+5) = 1.0;
      end

      colormap ( 'jet' );
      moviename = strcat ( 'exercise4_', num2str(f), '_', num2str(k), '.avi' );
      writer = VideoWriter ( moviename );
      open ( writer );

      figure( 1 );

    else

      uLaplace = laplacian9_torus ( u, dx, dy );
      vLaplace = laplacian9_torus ( v, dx, dy );

      dudt = ?
      dvdt = ?

      u = u + dt * dudt;
      v = v + dt * dvdt;

    end

    if ( mod ( it, frameRate ) == 1 )
      image ( u, 'CDataMapping', 'scaled' );
      label = sprintf ( 'Reaction-diffusion, step = %d\n', it );
      title ( label );
      axis ( 'equal' );
      axis ( 'tight' );
      set ( gca, 'xtick', [] );
      set ( gca, 'ytick', [] );
      frame = getframe ( 1 );
      writeVideo ( writer, frame );
    end

  end

  close ( writer );
  fprintf ( 1, '  Movie saved as "%s"\n', moviename );

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

  return
end
function L = laplacian9_torus ( A, dx, dy )

%*****************************************************************************80
%
%% laplacian9_torus uses a 9 point Laplacian stencil on a torus.
%
%  Discussion:
%
%    1) Assumes uniform, equal mesh in both directions.
%    2) Assumes the spacing is dx in both directions.
%
%  Modified:
%
%    28 June 2020
%
  [ nxp, nyp ] = size ( A );
  L = zeros ( nxp, nyp );

  for j = 1 : nyp

    for i = 1 : nxp

      L(i,j) = ?

    end
  end

  return
end

