Finite Difference Method for Poisson's Equation

Question

I'm writing a Python program to solve Poisson's equation

$$ \nabla^2 u = f \quad \mathrm{on} \quad \Omega $$ with $$ \frac{\partial u}{\partial n} = 0 \quad \mathrm{on} \quad \partial \Omega. $$

Here, I use a five-point stencil to discretize the Laplace operator and for each node in the computational domain, the resulting equation is

$$ \frac{1}{\delta x^2}(u_{i-1,j} + u_{i+1,j} - 2u_{i,j}) + \frac{1}{\delta y^2}(u_{i,j-1} + u_{i,j+1} - 2u_{i,j}) = f_{i,j}. $$

If there are $M$ cells in the $x$-direction and $N$-cells in the $y$-direction in the computational domain, I write my program so that each row of the corresponding matrix corresponds $only$ to the interior points of the grid - that is to say, $x(2:M-1)\times y(2:M-1)$.

import numpy as np
import matplotlib.pyplot as plt

plt.close('all')

# Limits in x
a = 0
b = 2

# Limits in y
c = 0
d = 1

# Computational cells in x-direction
M = 4

# Cells in y-direction
N=4
dx = (b - a)/M
dy=(d - c)/N
x=np.arange(a,2.0001,dx)
y=np.arange(c,1.0001,dy)
X, Y = np.meshgrid(x, y)

# Right-hand side function
def f(x,y):
    return np.cos(np.pi*x)*((12*y**2)-(12*y)+2-((np.pi**2)*(y**2)*(1-y)**2))

# Matrix only has entries for the (M-1)*(N-1) interior grid points 
# of the discretized computational domain
N2=(M-1)*(N-1)
A=np.zeros((N2,N2))

for j in range(N-1): 
    for i in range(M-1):
        modInd = int((j)*(M-1)) + i 
        # Handles corner nodes 
        if j == 0 and i == 0:
            A[modInd,modInd] = -(1/dx**2 + 1/dy**2)
        elif i == M-2 and j == 0:
            A[modInd,modInd] = -1/dx**2 - 1/dy**2
        elif i == 0 and j == N-2:
            A[modInd,modInd] = -1/dx**2 - 1/dy**2
        elif i == M - 2 and j == N-2:
            A[modInd,modInd] = -1/dx**2 - 1/dy**2 
            
        # Handles boundary nodes
        elif j == 0 and 0 < i < M -2:
            A[modInd,modInd] = -(2/dx**2 + 1/dy**2)
        elif j == N-2 and 0 < i < M - 2:
            A[modInd,modInd] = -(2/dx**2 + 1/dy**2)
        elif i == 0 and 0 < j < N - 2:
            A[modInd,modInd] = -(1/dx**2 + 2/dy**2)
        elif i == M-2 and 0 < j < N-2:
            A[modInd,modInd] = -(1/dx**2 + 2/dy**2)  
            
        # Interior nodes
        else:
            A[modInd][modInd] = -2*(1/dx**2 + 1/dy**2)
        
        #Handles off-diagonal matrix entries from
        #points left and right of node in 5-pt stencil
        if i >= 0 and i < M-2:
             A[modInd+1][modInd] = 1/dx**2  
        if i > 0 and i <= M-2:
             A[modInd-1][modInd] = 1/dx**2
            
        #Handles entries in identity sub-matrices from
        #points north and south of node in 5-pt stencil
        if j > 0 and j <=N-2:
                A[modInd][modInd - (M-1)] = 1/dy**2
        if j >=0 and j < N - 2:
                A[modInd][modInd + (M-1)] = 1/dy**2
 
r=np.zeros(N2)  
     
# Right-hand side vector r     
for i in range (0,M-1):
    for j in range (0,N-1):           
        r[i+(N-1)*j]=f(x[i+1],y[j+1])
        print(x[i+1],y[i+1])

r[0] = 0
A[0,:] = 0
A[0,0] = 1
        
U = np.linalg.solve(A,r)
U = U.reshape(N-1,N-1)

Xint = X[1:N,1:N]
Yint = Y[1:N,1:N]
[fig, ax] = plt.subplots()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(Xint,Yint,U)

The following should produce an approximation the is second order accurate, but it doesn't. Regardless of how fine I make the grid, my numerical solution is always downshifted by some amount. Any advice would be appreciated.