function [ U, X, Y, Z ] = fd3d_poisson ( nx, ny, nz, xmin, xmax, ymin, ymax, ...
  zmin, zmax, f, g )

%*****************************************************************************80
%
%% fd3d_poisson() solves a Poisson equation in a cube.
%
%  Discussion:
%
%    This program uses the finite difference method to estimate the solution
%    to a Poisson problem in a cube:
%
%      - Del U = f(x,y,z) in Omega
%            U = g(x,y,z) on dOmega
%
%    The user specifies 
%      the number of grid points (nx,ny,nz), 
%      the spatial extent (xmin,xmax,ymin,ymax,zmin,zmax),
%      the functions f(x,y,z), g(x,y,z).
%
%    No attempt is made to efficiently construct the matrix, or to store it sparsely.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    05 December 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer NX, NY, NZ, the number of grid points in the X, Y, and Z directions.
%
%    real XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX, the X, Y, and Z extents.
%
%    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,NZ), a table of solution values.
%
%    real X(NY,NX,NZ), Y(NY,NX,NZ), Z(NY,NX,NZ), the coordinates of grid points.
%

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

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

  return
end

