Eigenmodes of 2D Poisson matirix with Neumann boundary conditions

Question

I am trying to find analytic expressions for the eigenvectors (and eigenvalues) of the 2D discrete Poisson matrix, in the case of zero Neumann boundary conditions. [ pic 1 ]

In my case, I'm using a basic finite difference stencil for discretizing the 2D Poisson equation. The square mesh has N interior points in each direction (N = 3 in the figure). In this case there are n = N^2 = 9 mesh values $u_{ij}$ . For a node $(i, j)$, the stencil (and the boundary conditions) looks like:

$ 4u_{i,j} -u_{i,j+1} -u_{i,j-1}-u_{i+1,j}-u_{i-1,j} = \lambda u_{i,j} \quad $( for $i,j=2:N-1$)

Boundary conditions

$ u_{i,j} -u_{i,j+1}= \lambda u_{i,j} \quad (i=2:N-1,j=1)\\$
$ u_{i,j} -u_{i,j-1}= \lambda u_{i,j} \quad (i=2:N-1,j=N)\\$
$ u_{i,j} -u_{i+1,j}= \lambda u_{i,j} \quad (i=1,j=2:N-1)\\$
$ u_{i,j} -u_{i-1,j}= \lambda u_{i,j} \quad (i=N,j=2:N-1)\\$
$ u_{1,1} =u_{1,N}=u_{N,1} =u_{N,N} = 0\\$

Here, $u=u^{k,l}$ is the $(k,l)^{th}$ eigenvector with components $u_{i,j}$, where $1\le i,j,k,l \le n$.

[ pic 1 : 2D discrete Poisson matrix for $N*N = 9$]

[[ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  1.  0.  0. -1.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  1. -1.  0.  0.  0.  0.]
 [ 0. -1.  0. -1.  4. -1.  0. -1.  0.]
 [ 0.  0.  0.  0. -1.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0. -1.  0.  0.  1.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]]

I tried to figure out the analytic eigenmodes using discrete cosine transform which works for the 1D case (see eq 2. More precisely, what Wikipedia calls "DCT-2" works in the 1D Neumann case.)

[ DCT-2 ] Centers j = $-\frac{1}{2}$ and $N-\frac{1}{2}$ Components : $cos ((j+\frac{1}{2}k)\frac{\pi}{N})$

\begin{bmatrix} 1 & -1 & & &\\ -1& 2 &-1& &\\ & ...&& &\\ &-1&2&-1&\\ &&-1&1\\ \end{bmatrix}

However, that approach doesn't seem to work anymore in 2D. We tried the following choices: (eq 3 , eq 4)

The index $(i,j)$ represents the $i*dim+j$(flat index) eigenvalue or eigenvector and each runs from 1 to N. And index $k$ and $l$ represents component of an eigenvectors and rach runs from 1 to N.

$ U^{kl}_{ij} : cos((i-\frac{1}{2})\frac{k \pi}{N}))cos((j-\frac{1}{2})\frac{l \pi}{N})) \quad$ eq(3)

$ U^{kl}_{ij} : cos(i\frac{k \pi}{N}))cos(j\frac{l \pi}{N})) \quad $eq(4)

Is there any way to get analytic eigenvectors for this discrete Poisson system?