function bvp_fd ( n )

%*****************************************************************************80
%
%% bvp_fd() uses the finite difference method for a boundary value problem.
%
%  Discussion:
%
%    Demonstrate how the solution of a linear boundary value problem (BVP)
%    can be approximated using the finite difference method.
%
%    Choose n nodes in the interval [a,b].
%
%    Define n x n linear system.
%
%    First and last equations enforce boundary conditions.
%    Equations 2 through n-1 are discretized versions of the differential
%    equation.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    12 November 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer N: the number of nodes to use in the interval.
%

%
%  The BVP is specified over [a,b]:
%
  a = 0.0;
  b = 0.5;
%
%  The side conditions are y(a) = ya, y(b) = yb.
%
  ya = 1.0;
  yb = 2.0 * exp ( 0.5 );
%
%  Set the nodes.
%
  x = linspace ( a, b, n );
  dx = ( b - a ) / ( n - 1 );
%
%  Allocate the matrix and right hand side.
%
  A = zeros ( n, n );
  rhs = zeros ( n, 1 );
%
%  First equation implements a boundary condition.
%
  A(1,1) = 1.0;
  rhs(1) = ya;
%
%  Internal equations:
%  y" - y = 4 exp(x)
%
  for i = 2 : n - 1
    A(i,i-1) =  1.0 / dx^2; 
    A(i,i  ) = -2.0 / dx^2 - 1.0;
    A(i,i+1) =  1.0 / dx^2;
  end

  rhs(2:n-1) = 4.0 * exp ( x(2:n-1) );
%
%  Last equation implements a boundary condition.
%
  A(n,n) = 1.0;
  rhs(n) = yb;
%
%  Solve for y.
%
  y = A \ rhs;

  y_exact = ( 2.0 * x + 1.0 ) .* exp ( x );
%
%  Now plot the solution.
%
  clf ( );
  hold ( 'on' );
    plot ( x, y_exact, 'b-', 'linewidth', 2 );
    plot ( x, y, 'ro', 'linewidth', 2 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y(X) --->' );
  title ( 'Solution of y"=y+4e^x, y(0)=1, y(1/2)=2 sqrt(x).' );
  legend ( 'exact', 'FD solution' );
  filename = 'bvp_fd_test.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end

