You can use (for whatever dimension you like) where for 2D, h = {hx,hy}, n = {nx,ny}
zeroBoundaryLaplaceMatrix[h_List, n_List] /; Length[h] == Length[n] :=
Module[{d = Length[h], fdd, coords},
coords = MapThread[(#1 Range[0, #2]) &, {h, n}];
Total[Table[
fdd = NDSolve`FiniteDifferenceDerivative[2 UnitVector[d, i], coords, "DifferenceOrder" -> 2, "PeriodicInterpolation" -> True];
fdd["DifferentiationMatrix"], {i, d}]]]
The zero boundary points in each dimension are at 0 and h*(n + 1)
