function hw4 ( theta )

%*****************************************************************************80
%
%% hw4: theta method for a heat equation, with symmetry condition.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 July 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'hw4:\n' );
  fprintf ( 1, '  Use theta method for 1D heat equation\n' );
  fprintf ( 1, '  with symmetry boundary condition.\n' );
%
%  Space
%
  nx = 21;
  a = 0.0;
  b = 1.0;
  x = linspace ( a, b, nx );
  dx = ( b - a ) / ( nx - 1 );
%
%  Time
%
  nt = 11;
  c = 0.0;
  d = 1.0;
  t = linspace ( c, d, nt );
  dt = ( d - c ) / ( nt - 1 );
%
%  Heat coefficient.
%
  k = 4.0 / pi^2;
%
%  Compute the CFL parameter.
%  If theta < 1/2, then we can have problems if cfl > 1/2!
%
  cfl = k * dt / dx^2;
  fprintf ( 1, '  CFL parameter is %g\n', cfl );
%
%  Create solution array, set initial condition.
%
  U = zeros ( nx, nt );
  U(1:nx,1) = g ( x, t(1) );
  umin = 0.0;
  umax = 1.0;
%
%  Create matrix and right hand side.
%
  A = zeros ( nx, nx );
  rhs = zeros ( nx, 1 );
%
%  Plot solution j, set up and solve system for solution j+1.
%
  for j = 1 : nt - 1
%
%  Show snapshot of solution and pause.
%
    plot ( x, U(1:nx,j), 'ro', x, g(x,t(j) ), 'b-', 'linewidth', 3 );
    axis ( [ a, b, umin, umax ] );
    grid ( 'on' );
    label = sprintf ( '%g method+Symmetry, Time step %d', theta, j );
    title ( label );
    pause ( )
%
%  Update solution.
%
    A(1,1) = 1.0;
    rhs(1) = g ( x(1),t(j+1) );

    for i = 2 : nx - 1
      A(i,i-1) =     - theta * 1.0 * k * dt / dx^2;
      A(i,i)   = 1.0 + theta * 2.0 * k * dt / dx^2;
      A(i,i+1) =     - theta * 1.0 * k * dt / dx^2;
      rhs(i) = U(i,j) + ( 1.0 - theta ) * k * dt / dx^2 ...
        * ( U(i-1,j) - 2.0 * U(i,j) + U(i+1,j) );
    end
%
%  Symmetry condition at right endpoint allows us to "reflect" coefficient
%  of U(nx+1,j) to U(nx-1,j), and to replace U(i+1,j) by U(-1,j) in right hand side.
%
    i = nx;
    A(i,i-1) =     - theta * 2.0 * k * dt / dx^2;
    A(i,i)   = 1.0 + theta * 2.0 * k * dt / dx^2;
    rhs(i) =      U(i,j) + ( 1.0 - theta ) * k * dt / dx^2 ...
      * ( U(i-1,j) - 2.0 * U(i,j) + U(i-1,j) );

    U(1:nx,j+1) = A \ rhs;

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

  return
end
function value = g ( x, t )

%*****************************************************************************80
%
%% g evaluates the exact solution.
%
% x = x(:);
  value = sin ( pi * x / 2.0 ) .* exp ( - t );

  return
end

