function test03 ( )

%*****************************************************************************80
%
%% test03 uses a discontinuous initial condition.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    29 March 2021
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test03:\n' );
  fprintf ( 1, '  Compute an approximate solution to the time-dependent\n' );
  fprintf ( 1, '  one dimensional heat equation:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    dH/dt - K * d2H/dx2 = 0\n' );
  fprintf ( 1, '\n' );
%
%  Heat coefficient.
%
  k = 2.0;
  fprintf ( 1, '  Constant K = %g\n', k );
%
%  Compute equally spaced spatial nodes X.
%
  x_num = 21;
  x_min = -5.0;
  x_max = +5.0;
  dx = ( x_max - x_min ) / ( x_num - 1 );
  x = linspace ( x_min, x_max, x_num );
  fprintf ( 1, '  Number of X nodes = %d\n', x_num );
  fprintf ( 1, '  X interval is [%f,%f]\n', x_min, x_max );
  fprintf ( 1, '  X spacing is %f\n', dx );
%
%  Compute equally spaced time nodes T.
%
  t_num = 81;
  t_min = 0.0;
  t_max = 2.0;
  dt = ( t_max - t_min ) / ( t_num - 1 );
  t = linspace ( t_min, t_max, t_num );
  fprintf ( 1, '  Number of T values = %d\n', t_num );
  fprintf ( 1, '  T interval is [%f,%f]\n', t_min, t_max );
  fprintf ( 1, '  T spacing is %f\n', dt );
%
%  Get the CFL coefficient.
%
  cfl = k * dt / dx^2;
  fprintf ( 1, '  CFL coefficient = %g\n', cfl );
%
%  Get the system matrix.
%  First and last equation are Dirichlet boundary conditions.
%  Intermediate equations are a version of:
%    H(new) - H(old)       H(new,left) - 2 H(new) + H(new,right)
%    --------------- - K * ------------------------------------- = 0
%          dt                              dx^2
%  after multiplying by dt and moving H(old) to the right hand side.
%
  a = sparse ( [], [], [], x_num, x_num );

  a(1,1) = 1.0;

  for i = 2 : x_num - 1
    a(i,i-1) =           - cfl;
    a(i,i  ) = 1.0 + 2.0 * cfl;
    a(i,i+1) =           - cfl;
  end

  a(x_num,x_num) = 1.0;
%
%  Repeatedly solve the linear system.
%
  h = zeros ( x_num, 1 );
  hmat = zeros ( x_num, t_num );

  for j = 1 : t_num
%
%  H(new) = initial condition.
%
    if ( j == 1 )

      i1 = ( x(1:x_num) <=  5.0 );
      i2 = ( x(1:x_num) <   1.0 );
      i3 = ( x(1:x_num) <  -1.0 );
      h(i1,1) = 25.0;
      h(i2,1) =  0.0;
      h(i3,1) = 15.0;
%
%  Solve A * H(new) = H(old)
%
    else
      h = a \ h;
    end

    hmat(1:x_num,j) = h(1:x_num);

  end
%
%  Plot the data.
%
  figure ( 3 )
  [ xmat, tmat ] = ndgrid ( x, t );
  mesh ( xmat, tmat, hmat );
  title ( 'H(X,T): test03' );
  xlabel ( '<-- X -->' );
  ylabel ( '<-- T -->' );
  zlabel ( '<--  H(X,T) -->' );
  filename = 'test03.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end
