function exercise4 ( )

%*****************************************************************************80
%
%% exercise4 solves the time-dependent heat equation with insulated ends.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 May 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise4\n' );
  fprintf ( 1, '  Time dependent heat equation with insulation.\n' );
%
%  Set the X values.
%
  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 the first solution to the initial condition.
%
  U = zeros ( nx, nt );

  U(1:nx,1) = h ( x );
  U(1,1) = U(2,1);
  U(nx,1) = U(nx,2);

  E = dx * ( sum ( U(1:nx,1) ) - 0.5 * U(1,1) - 0.5 * U(nx,1) )
%
%  Compute values at time j+1 from values at time j.
%  Enforce the Neumann boundary conditions.
%  Estimate the energy E.
%
  for j = 1 : nt - 1
%
%  Show snapshot of solution and pause.
%
    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 ( )
%
%  Update solution.
%
    U(2:nx-1,j+1) = U(2:nx-1,j) ...
      + ( dt / dx^2 ) * ( U(1:nx-2,j) - 2.0 * U(2:nx-1,j) + U(3:nx,j) );

    U(1,j+1) = U(2,j+1);
    U(nx,j+1) = U(nx-1,j+1);

    E = dx * ( sum ( U(1:nx,j+1) ) - 0.5 * U(1,j+1) - 0.5 * U(nx,j+1) ); 

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise4\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

