function l = laplacian5_plane ( u, dx, dy )

%*****************************************************************************80
%
%% laplacian5_plane() approximates the laplacian on a plane grid.
%
%  Discussion:
%
%    The domain is represented by a grid of nx * ny points.
%    There is constant spacing dx and dy in the x and y directions respectifely.
%
%    The laplacian is computed for the interior points, but a value
%    of 0 is returned for points on the boundary.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    25 September 2024
%
%  Input:
%
%    real u(nx,ny): the function values at grid points.
%
%    real dx, dy: the spacing in the x and y directions.
%
%  Output:
%
%    real l(nx,ny): the estimate for the laplacian.
%
  [ nx, ny ] = size ( u );

  uxx = zeros ( nx, ny );
  uyy = zeros ( nx, ny );

  uxx(2:nx-1,:) = ( u(3:nx,:) - 2.0 * u(2:nx-1,:) + u(1:nx-2,:) ) / dx^2;
  uyy(:,2:ny-1) = ( u(:,3:ny) - 2.0 * u(:,2:ny-1) + u(:,1:ny-2) ) / dy^2;

  l = uxx + uyy;

  return
end

