function exercise5 ( )

%*****************************************************************************80
%
%% exercise5 solves the time-dependent heat equation with periodic conditions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 May 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise4\n' );
  fprintf ( 1, '  Time dependent heat equation with periodic boundary.\n' );
%
%  Set the X values.
%  The last node is the same as the first node.
%  So we really only have NX-1 nodes!
%
  nx = 21;
  a = 0.0;
  b = 3.0;
  x = linspace ( a, b, nx );
  dx = ( b - a ) / ( nx - 1 );
%
%  Set the T values.
%
  nt = 101;
  c = 0.0;
  d = 1.0;
  t = linspace ( c, d, nt );
  dt = ( d - c ) / ( nt - 1 );

  fprintf ( 1, '  Use nx = %d and nt = %d.\n', nx, nt );
%
%  Set space for the solution.
%  The value at node NX is always a copy of the value at node 1.
%
  U = zeros ( nx, nt );
%
%  Make it easy to refer to left, center, and right nodes.
%
  l = [ nx - 1, 1 : nx - 2 ];
  c = 1 : nx - 1;
  r = [ 2 : nx-1, 1 ];
%
%  Set the first solution to the initial condition.
%
  j = 1;
  U(c,j) = h ( x(c) );
  U(nx,j) = U(1,j);
  E = dx * ( sum ( U(c,j) ) );

  for j = 1 : nt - 1
%
%  Plot solution j
%
    plot ( x, U(1:nx,j) )
    axis ( [ 0.0, 3.0, 0.0, 3.0 ] )
    grid ( 'on' );
    label = sprintf ( 'Time step %d, E = %g', j, E );
    title ( label );
    pause ( )
%
%  Compute solution j+1.
%
    U(c,j+1) = U(c,j) + ( dt / dx^2 ) * ( U(l,j) - 2.0 * U(c,j) + U(r,j) );
    U(nx,j+1) = U(1,j+1);

    E = dx * ( sum ( U(c,j+1) ) ); 

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

  return
end
function value = h ( x )

%*****************************************************************************80
%
%% H evaluates the initial condition.
%
  value =   x .* ( 3 - x ) .* ( x.^2 - 2*x + 1 );

  return
end

