function [ U, X, Y ] = fd2d_poisson ( nx, ny, xmin, xmax, ymin, ymax, f, g )

%*****************************************************************************80
%
%% fd2d_poisson is a simple program to solve a Poisson equation in a rectangle.
%
%  Discussion:
%
%    This program uses the finite difference method to estimate the solution
%    to a Poisson problem on a rectangle:
%
%      - Del U = f(x,y) in Omega
%            U = g(x,y) on dOmega
%
%    The user specifies the number of grid points (nx,ny), the spatial
%    extent (xmin,xmax,ymin,ymax), and the functions f(x,y), g(x,y).
%
%    No attempt is made to efficiently construct the matrix, or to store
%    it sparsely.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    06 December 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NX, NY, the number of grid points in the X and Y directions.
%
%    real XMIN, XMAX, YMIN, YMAX, the X and Y endpoints.
%
%    function handle F, the name of the MATLAB function which evaluates
%    the right hand side of the Poisson equation.
%
%    function handle G, the name of the MATLAB function which evaluates
%    the Dirichlet boundary conditions.
%
%  Output:
%
%    real U(NY,NX), a table of solution values.
%
%    real X(NY,NX), Y(NY,NX), tables of the X and Y coordinates of the points.
%

%
%  Compute the X and Y spacings.
%
  hx = ( xmax - xmin ) / ( nx - 1 );
  hy = ( ymax - ymin ) / ( ny - 1 );
%
%  Compute the X and Y coordinate vectors.
%
  x = linspace ( xmin, xmax, nx );
  y = linspace ( ymin, ymax, ny );
%
%  Set the system matrix A and right hand side rhs.
%
  A = zeros ( nx * ny, nx * ny );
  rhs = zeros ( nx * ny, 1 );

  eqn = 0;
  for j = 1 : ny
    for i = 1 : nx
      eqn = eqn + 1;
      if ( i == 1 || i == nx || j == 1 || j == ny )
        A(eqn,eqn) = 1.0;
        rhs(eqn) = g ( x(i), y(j) );
      else
        A(eqn,eqn-nx) = - 1.0 / hy^2;
        A(eqn,eqn-1)  = - 1.0 / hx^2;
        A(eqn,eqn)    =   2.0 / hx^2 + 2.0 / hy^2;
        A(eqn,eqn+1)  = - 1.0 / hx^2;
        A(eqn,eqn+nx) = - 1.0 / hy^2;
        rhs(eqn)      = f ( x(i), y(j) );
      end
    end
  end
%
%  Solve for u, the solution vector.
%
  u = A \ rhs;
%
%  Return matrices U, X, Y.
%
  U = reshape ( u, nx, ny );
  [ X, Y ] = ndgrid ( x, y );

  return
end

