function [ X, Y, U ] = poisson_solver ( a, b, c, d, nx, ny, f, g )

%*****************************************************************************80
%
%% poisson_solver solves a 2D Poisson equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 July 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real a, b, c, d: the x and y ranges.
%
%    integer nx, ny: the number of points in the x and y directions.
%
%    function value = f(x,y): the right hand side.
%
%    function value = g(x,y): the exact solution.
%
%  Output:
%
%    real X(nx,ny), Y(nx,ny): tables of coordinate values.
%
%    real U(nx,ny): a table of the solution values.
%
  hx = ( b - a ) / ( nx - 1 );
  x = linspace ( a, b, nx );
  hy = ( d - c ) / ( ny - 1 );
  y = linspace ( c, d, ny );

  A = zeros(nx*ny,nx*ny);  
  rhs = zeros(nx*ny,1);

  k = 0;
  for j = 1 : ny
    for i = 1 : nx
      k = k + 1;
      if ( j == 1 || j == ny || i == 1 || i == nx )
        A(k,k) = 1.0;
        rhs(k) = g ( x(i), y(j) );
      else
        CENTER = k;
        A(k,CENTER) = 2.0 / hx^2 + 2.0 / hy^2;
        WEST = k_from_ij ( i-1, j, nx );
        A(k,WEST) = -1.0 / hx^2;
        EAST = k_from_ij ( i+1, j, nx );
        A(k,EAST) = -1.0 / hx^2;
        SOUTH = k_from_ij ( i, j-1, nx );
        A(k,SOUTH) = -1.0 / hy^2;
        NORTH = k_from_ij ( i, j+1, nx );
        A(k,NORTH) = -1.0 / hy^2;
        rhs(k) = f ( x(i), y(j) );
      end
    end
  end

  u = A \ rhs;
%
%  Create tables for output.
%
  U = reshape ( u, nx, ny );
  [ X, Y ] = ndgrid ( x, y );

  return
end

