function exercise3 ( )

%*****************************************************************************80
%
%% exercise3: Crank Nicolson method for a heat equation with Neumann condition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 July 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise3:\n' );
  fprintf ( 1, '  Use Crank-Nicolson method for 1D heat equation\n' );
  fprintf ( 1, '  with Neumann boundary condition.\n' );
%
%  Space
%
  nx = 21;
  a = 0.0;
  b = 1.0;
  x = linspace ( a, b, nx );
  dx = ( b - a ) / ( nx - 1 );
%
%  Time
%
  nt = 11;
  c = 0.0;
  d = 1.0;
  t = linspace ( c, d, nt );
  dt = ( d - c ) / ( nt - 1 );
%
%  Heat coefficient.
%
  k = 4.0 / pi^2;
%
%  Compute the CFL parameter.
%  Because we are using Crank Nicolson, we don't care if cfl > 1/2!
%
  cfl = k * dt / dx^2;
  fprintf ( 1, '  CFL parameter is %g\n', cfl );
%
%  Create solution array, set initial condition.
%
  U = zeros ( nx, nt );
  U(1:nx,1) = g ( x, t(1) );
  umin = 0.0;
  umax = 1.0;
%
%  Create matrix and right hand side.
%
  A = zeros ( nx, nx );
  rhs = zeros ( nx, 1 );
%
%  Plot solution j, set up and solve system for solution j+1.
%
  for j = 1 : nt - 1
%
%  Show snapshot of solution and pause.
%
    plot ( x, U(1:nx,j), 'ro', x, g(x,t(j) ), 'b-', 'linewidth', 3 );
    axis ( [ a, b, umin, umax ] );
    grid ( 'on' );
    label = sprintf ( 'Crank Nicolson/Neumann, Time step %d', j );
    title ( label );
    pause ( )
%
%  Update solution.
%
    A(1,1) = 1.0;
    rhs(1) = g ( x(1),t(j+1) );

    for i = 2 : nx - 1
      A(i,i-1) =     - 1.0 * k * dt / 2.0 / dx^2;
      A(i,i)   = 1.0 + 2.0 * k * dt / 2.0 / dx^2;
      A(i,i+1) =     - 1.0 * k * dt / 2.0 / dx^2;
      rhs(i) = U(i,j) + k * dt / 2.0 / dx^2 * ( U(i-1,j) - 2.0 * U(i,j) + U(i+1,j) );
    end
%
%  Neumann boundary condition at right.
%
    i = nx;
    A(i,i-1) = -1.0 / dx;
    A(i,i)   = +1.0 / dx;
    rhs(i) = 0.0

    U(1:nx,j+1) = A \ rhs;

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise3\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function value = g ( x, t )

%*****************************************************************************80
%
%% g evaluates the exact solution.
%
% x = x(:);
  value = sin ( pi * x / 2.0 ) .* exp ( - t );

  return
end

