function D = delsq ( G )

%*****************************************************************************80
%
%% delsq() constructs the five-point finite difference Laplacian.
%
%  Discussion:
%
%    delsq(G) is the sparse form of the two-dimensional,
%    5-point discrete negative Laplacian on the grid G.
%
%  Licensing:
%
%    Based on code that is copyright 1984-2013 The MathWorks, Inc. 
%
%  Modified:
%
%    30 August 2023
%
%  Author:
%
%    Original MATLAB version by Cleve Moler.
%    This version by John Burkardt
%
%  Input:
%
%    integer g(m,m): a grid whose positive entries are the indices of
%    nodes.
%
%  Output:
%
%    sparse matrix D(): defines the negative Laplacian operator 
%    for the grid.
%
  m = size ( G, 1 );
%
%  Interior points have a nonzero value in G.
%
  p = find ( G );
%
%  Connect interior points to themselves with 4's.
%
  i = G(p);
  j = G(p);
  s = 4.0 * ones ( size ( p ) );
%
%  Using 1-dimensional array index, the neighbors of node p are
%  p-1 (north)
%  p+m (east)
%  p+1 (south)
%  p-m (west)
%  Each of these neighbors gets a coefficient of -1.
%
  for k = [ -1, m, 1, -m ]
    Q = G(p+k);
    q = find ( Q );
    i = [ i; G(p(q)) ];
    j = [ j; Q(q) ];
    s = [ s; -ones(length(q),1) ];
  end
%
%  Define D as a sparse matrix.
%
  D = sparse ( i, j, s );

  return
end

