function exercise1 ( )

%*****************************************************************************80
%
%% exercise1 solves the time-independent heat equation.
%
%  Discussion:
%
%    This program seeks a solution of the steady heat equation:
%
%      - d/dx dUdx = f(x)
%
%    over the interval [0,1] with 
%
%      f(x) = 100 for 0.5 <= x <= 0.7
%               0 otherwise
%
%    with Dirichlet boundary conditions
%
%      U(0) = 0.0
%      U(1) = 5.0.
%
%    The code uses the finite difference method to approximate the
%    second derivative in space.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise1\n' );
  fprintf ( 1, '  Steady heat equation with nonzero right hand side.\n' );

  n = 41;
  a = 0.0;
  b = 1.0;
%
%  Set the X values.
%
  x = linspace ( a, b, n );
  dx = ( b - a ) / ( n - 1 );

  for heating = [ false, true ]
%
%  Set up the matrix.
%
    A = zeros ( n, n );
%
%  Left Dirichlet boundary condition.
%
    A(1,1) = 1.0;
    rhs(1,1) = 0.0;

    for i = 2 : n - 1
      A(i,i-1) = - 1.0 / dx^2;
      A(i,i) =     2.0 / dx^2;
      A(i,i+1) = - 1.0 / dx^2;
      rhs(i,1) = f ( x(i), heating );
    end
%
%  Right Dirichlet boundary condition.
%
    A(n,n) = 1.0;
    rhs(n,1) = 80.0;
%
%  Solve the linear system.
%
    u = A \ rhs;
%
%  Plot the result.
%
    figure ( );
    hold ( 'on' );
    plot ( x, u, 'r.', 'markersize', 20 );
    umax = max ( u );
    plot ( [ 0.5, 0.5 ], [0.0, umax], 'k--' );
    plot ( [ 0.7, 0.7 ], [0.0, umax], 'k--' );
    hold ( 'off' );
    grid ( 'on' );
    xlabel ( '<-- X -->' );
    ylabel ( '<-- T(X) -->' );
    if ( heating )
      title ( 'exercise1 with heating' );
      filename = 'exercise1_hot.png';
    else
      title ( 'exercise1 with no heating' );
      filename = 'exercise1_cold.png';
    end
    print ( '-dpng', filename );
    fprintf ( 1, '  Graphics saved as "%s"\n', filename );

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

  return
end
function value = f ( x, heating )

%*****************************************************************************80
%
%% F evaluates the right hand side function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, the position.
%
%    logical HEATING: true if an internal heat source is applied.
%
%  Output:
%
%    real VALUE, the right hand side function.
%
  if ( heating )
    value = 1500.0 * ( 0.5 <= x && x <= 0.7 );
  else
    value = 0.0;
  end

  return
end

