function hw3 ( )

%*****************************************************************************80
%
%% hw3 solves a version of the heat equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2020
%
%  Author:
%
%    John Burkardt
%
  nx = 21;
  nt = 11;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'hw3\n' );
  fprintf ( 1, '  Time dependent heat equation.\n' );
  fprintf ( 1, '  nx = %d, nt = %d\n', nx, nt );
%
%  Set the X values.
%
  a = 0.0;
  b = pi;
  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 );
%
%  Set the heat coefficient k.
%
  k = 1.0 / 9.0;
%
%  Report the CFL parameter.
%
  cfl = k * 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, t(1) );
%
%  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) );
%
%  Define L, C, R for left, center, and right sets of nodes.
%
  l = 1 : nx-2;
  c = 2 : nx-1;
  r = 3 : nx;
%
%  Use values at time J to compute values at time j+1.
%
  for j = 1 : nt - 1
%
%  Show snapshot of solution and pause.
%
    V = h ( x, t(j) );
    plot ( x, U(1:nx,j), 'ro', x, V, 'b-', 'linewidth', 3 );
    axis ( [ a, b, umin, umax ] )
    grid ( 'on' );
    label = sprintf ( 'Time step %d, t = %g', j, t(j) );
    title ( label );
    pause ( )
%
%  Update solution.
%
    U(1,j+1) = h ( x(1), t(j+1) );

    U(c,j+1) = U(c,j) + ( k * dt / dx^2 ) * ( U(l,j) - 2.0 * U(c,j) + U(r,j) );

    U(nx,j+1) = h ( x(nx), t(j+1) );

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

  return
end
function value = h ( x, t )

%*****************************************************************************80
%
%% H evaluates the exact solution.
%
%  Discussion:
%
%    dHdt   = - 3.0 * sin ( 3.0 * x ) * exp ( -t ) 
%    d2Hdx2 = -27.0 * sin ( 3.0 * x ) * exp ( -t )
%
%    Heat equation is dHdt = 9 * d2Hdx2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real x, t: the space and time coordinates.
%
%  Output:
%
%    real value: the value of the exact solution at x and t.
%
  value = 3.0 * sin ( 3.0 * x ) * exp ( - t );

  return
end


