FFT based Fast Poisson Solver convergence error rate giving an order below of expected

Question

Hello to the stackexchange community. I have been assigned a homework regarding an FFT based Poisson solver. However, I have been struggling with the error of convergence of this method. Before I ask my question I would like to give some context to the problem.

Take the Poisson equation

$$-\nabla^2 u = f$$

which we discretize as

$$-\frac{u_{i,j+1} -2u_{i,j} + u_{i,j-1}}{h^2} - \frac{u_{i+1,j} -2u_{i,j} + u_{i-1,j}}{h^2} = f_{i,j}.$$

We can rewrite the above in the form

$$AU + UA = F,$$

where $A= \text{diag}(-1,2,-1)$ is a symmetric tridiagonal matrix. We can decompose the matrix as $A=V\Lambda V^{-1}$, so

\begin{align} (V\Lambda V^{-1})U + U(V\Lambda V^{-1}) &= F \\ \implies \Lambda (V^{-1}UV) + (V^{-1}UV)\Lambda &= V^{-1}FV \\ \implies \Lambda \bar{U} + \bar{U}\Lambda &= \bar{F} \end{align}

where $\bar{U} = V^{-1}UV$, $\bar{F} = V^{-1}FV$, $\Lambda = \text{diag}(\lambda_i)$ is the eigenvalue matrix with $\lambda_i = -4\sin^2\left(\frac{i\pi}{2(N+1)}\right)$, and $V$ is the eigenvector matrix with elements $v_{i,j} = \sqrt{\frac{2}{N+1}}\sin\left(\frac{ij\pi}{N+1}\right)$ with $V^{-1} = \frac{2}{N+1}V$.

Now, the matrix multiplication $FV$ is the same as doing a discrete sine transform (DST) into the rows of $F$, while $V^{-1}F$ (I suppose) is the same as applying DST into the columns of $F$.

To get $U$ back, we perform

$$\bar{u}_{i,j} = \frac{\bar{f}_{i,j}}{\lambda_i+\lambda_j}$$

and

$$U = V\bar{U}V^{-1}$$

With that, we have our fast Poisson Solver. To get the rate of convergence of the method, we define it as

$$c = \frac{\log(\frac{e_1}{e_2})}{\log{2}}$$

where $e_1$, and $e_2$ are the errors with spacing $h/2$ and $h$, respectively. If $c$ is close to $2$, it would then indicate that the method is of order/rate $2$.

So, I define my problem as

$$f = -2\pi^2\sin(\pi x)\sin(\pi y),$$ for $0\leq x,y \leq 1$ with boundary conditions equal to $0$ at all boundaries. When I try to solve the problem, I get a convergence value of 1 instead of 2.

Therefore, I do not know what is going on. Below I attached a simple code I created for problem solution (some of the constants multiplying after the DST are a normalization).

import numpy as np
import scipy.fftpack as fft
import matplotlib.pyplot as plt

def anaf(x,y):
  varo = np.sin(np.pi*x)*np.sin(np.pi*y);
  return varo

inds = np.arange(5,11+1);
errorm = np.zeros(inds.shape)
for i,ind in enumerate(inds):
  n = 2**ind;
  x = np.linspace(0,1,n);
  y = np.linspace(0,1,n);
  h = x[1]-x[0];

  xg,yg = np.meshgrid(x,y);

  uana = anaf(xg,yg)
  F = -2*np.pi**2*uana;

  ij = np.arange(1,n+1);
  ik = np.arange(1,n+1);

  ijg,ikg = np.meshgrid(ij,ik);

  V = np.matrix(np.sqrt(2/(n+1))*(np.sin(ijg*ikg*np.pi/(n+1))));

  D  = -4/h**2 * (np.sin((ij*np.pi) / (2*(n + 1))))**2;

  F2 = .5*fft.dst(F.T,n=n,type=1);
  F2 = .5*fft.dst(F2.T,n=n,type=1);
  unum = np.zeros([n,n]);
  for j in range(len(x)):
    for k in range(len(y)):
      unum[j,k] = (2/n+1)**2*F2[j,k]/(D[j]+D[k]);

  unum = .5*fft.dst(unum.T,n=n,type=1);
  unum = .5*(2/(n+1))**2*fft.dst(unum.T,n=n,type=1);
  errorm[i] = np.abs(unum-uana).max();
  print('%s: %1.4f' % (ind,errorm[i]))

print('Conv: '+str(np.log(errorm[:-1]/errorm[1:])/np.log(2)));

plt.scatter(inds,np.log2(errorm));
plt.grid();
plt.show();

Please, excuse me if I made any mistakes or if made the question a bit confusing.