Discretisation of 2D Poisson equation in Sylvester matrix form

I am working on fast Poisson solvers and I have to understand some basic concepts with the discretisation of the 2D Poisson equation ($$U_{xx}+U_{yy}=f$$) in the Sylvester equation form ($$KX+XK=F$$).

I have been able to discretise the 2D Poisson equation using the central difference approximation (assuming $$\Delta x=\Delta y=h$$) and I got $$-\dfrac{1}{h^{2}}\left[ 2U_{i,j}-U_{i-1,j}-U_{i+1,j}\right]-\dfrac{1}{h^{2}}\left[ 2U_{i,j}-U_{i,j-1}-U_{i,j+1}\right]=f_{i,j}$$ Now, instead of simplifying this to get a linear equation $$A\vec{u}=\vec{f}$$ where $$A$$ is an $$(n-1)(n-1)\times (n-1)(n-1)$$ matrix and $$\vec{u}, \vec{f}$$ are column vectors. I need to express the discretised 2D Poisson equation in the Sylvester equation form $$KX+XK=F$$ where $$K=-\dfrac{1}{h^{2}}\begin{bmatrix} 2&-1&0&0&0&\ldots&0\\ -1&2&-1&0&0&\ldots&0\\ 0&-1&2&-1&0&\ldots&0\\ \vdots&\ddots&\ddots&\ddots&\ddots&\ddots&\vdots\\ 0&\ldots&0&-1&2&-1&0\\ 0&\ldots&\ldots&0&-1&2&-1\\ 0&\ldots&\ldots&\ldots&0&-1&2\\ \end{bmatrix}$$ and matrix $$X$$ represents the values of the solution on the interior nodes of a $$(n+1)\times (n+1)$$ equispaced grid.

I finding it difficult to find the transition from the discretised 2D Poisson equation to the Sylvester equation. Please I need help from whoever has idea about this.