 Discrete sine transform Poisson solver

For $-\Delta u=f$ on a square with homogeneous Dirichlet boundary conditions, the five-point finite difference operator on a mesh $h=1/N$ is diagonalized by $\sin (\pi ki/N)\sin (\pi \ell j/N)$. Its eigenvalues are ${\lambda }_{k\ell }=4h^{-2}[{\sin }^2(\pi k/(2N))+{\sin }^2(\pi \ell /(2N))]$. A two-dimensional discrete sine transform, division by these nonzero eigenvalues and the inverse transform give the solution. Odd extension realizes each one-dimensional transform by a Fast Fourier transform of length $2N$, making the total cost $O(N^2\log N)$ and storage $O(N^2)$. These finite-difference eigenvalues should not be confused with the continuous spectral eigenvalues ${\pi }^2(k^2+{\ell }^2)$.