function exercise3 ( nx, nt )

%*****************************************************************************80
%
%% exercise3 accepts nx and nt as parameters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 May 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NX, NT: the number of grid points in space and in time.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise3\n' );
  fprintf ( 1, '  Time dependent heat equation.\n' );
  fprintf ( 1, '  User specifies nx = %d, nt = %d\n', nx, nt );
%
%  Set the X values.
%
  a = 0.0;
  b = 3.0;
  x = linspace ( a, b, nx );
  dx = ( b - a ) / ( nx - 1 );
%
%  Set time values.
%
  c = 0.0;
  d = 1.0;
  t = linspace ( c, d, nt );
  dt = ( d - c ) / ( nt - 1 );
%
%  Report the CFL parameter.
%
  cfl = dt / dx^2;

  if ( 0.5 <= cfl )
    fprintf ( 1, '  The CFL parameter %g exceeds 1/2!\n', cfl );
    fprintf ( 1, '  The approximation method is unstable.\n' );
  else
    fprintf ( 1, '  The CFL parameter %g is below 1/2.\n', cfl );
    fprintf ( 1, '  The approximation method is stable.\n' );
  end
%
%  Set the first solution to the initial condition.
%
  U = zeros ( nx, nt );

  U(1:nx,1) = h ( x );
%
%  Record min and max of initial condition, so all plots use same scale.
%
  umin = min ( U(1:nx,1) );
  umax = max ( U(1:nx,1) );
%
%  Use values at time J to compute values at time j+1.
%
  for j = 1 : nt - 1
%
%  Show snapshot of solution and pause.
%
    plot ( x, U(1:nx,j) )
    axis ( [ a, b, umin, umax ] )
    grid ( 'on' );
    label = sprintf ( 'Time step %d', j );
    title ( label );
    pause ( )
%
%  Update solution.
%
    U(1,j+1) = h ( x(1) );
    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(nx,j+1) = h ( x(nx) );

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

