function bvp_fd_test ( )

%*****************************************************************************80
%
%% bvp_fd_test uses the finite difference method for a boundary value problem.
%
%  Discussion:
%
%    Demonstrate how the solution of a lineary 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 ODE.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    18 August 2019
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'bvp_fd_test\n' );
  fprintf ( 1, '  A simple example of the finite difference method,\n' )
  fprintf ( 1, '  in which a boundary value problem (BVP) is solved\n' );
  fprintf ( 1, '  by choosing a grid of n points, setting up a linear system\n' );
  fprintf ( 1, '  which enforces boundary conditions, and uses discretized\n' );
  fprintf ( 1, '  versions of derivatives at each point.\n' );
%
%  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.
%
  n = 11;
  x = linspace ( a, b, n );
  h = ( b - a ) / ( n - 1 );
%
%  Define the matrix and right hand side
%
  A = zeros ( n, n );
  rhs = zeros ( n, 1 );
%
%  First equation.
%
  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 / h^2; 
    A(i,i  ) = -2.0 / h^2 - 1.0;
    A(i,i+1) =  1.0 / h^2;
  end

  rhs(2:n-1) = 4.0 * exp ( x(2:n-1) );
%
%  Last equation.
%
  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.
%
  plot ( x, y_exact, 'b-', x, y, 'ro' );
  grid ( 'on' );
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y(X) --->' );
  title ( 'Solution of y"=y+4e^x, y(0)=1, y(1/2)=2 sqrt(x).' )
  filename = 'bvp_fd_test.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'bvp_fd_test\n' );
  fprintf ( 1, '  Normal end of execution.\n' )
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

