Implementing Finite Differences solver for 2D Poisson Equation

Question

I'm trying to solve the 2D Poisson equation:

$$ \begin{cases} -\Delta u = f & \text{in} \hspace{0.2cm} \Omega=(0,1)^{2} \\ g = u & \text{on} \hspace{0.2cm}\partial \Omega \end{cases} $$

Using finite differences with order of accuracy 2, by discretizing the square into a mesh with step size $h$ in each direction. I chose $u$ and $f$ such that:

\begin{align*} u(x,y) &= sin(\pi x) \\ f(x,y) &= \pi^2 sin(\pi x) \end{align*}

I saw that this could be written as a linear system $A_{h}^{2D} \mathbf{u} = \mathbf{f}$ where $\mathbf{u}$ is the vectorization of the matrix $(U_{h})_{i,j} = (u_{h}(x_{i}, y_{j}))$ and $\mathbf{f}$ is the vectorization of $(F)_{i,j} = (f(x_{i}, y_{j}))$, and:

$$ A_{h}^{2D} = I \otimes A_{h} + A_{h} \otimes I $$

where $A_{h} = \frac{1}{h^{2}}tridiag(-1,2,-1)$ and $\otimes$ denotes the Kronecker product. The following is my attempt at an implementation, however the actual solution (right) is very far off from the computed one (left) even as I make the step size h very small. I'm wondering if anyone can tell what could be going wrong in the implementation in Octave (MATLAB)? Thanks so much.

clear all
close all

f = @(x,y) (pi^2)*sin(pi*x);
u = @(x,y) sin(pi*x);

N = 20;
h = 1/(N+1);

Ah = (1/(h^2))*(diag(-1*ones(N-1,1),-1) + diag(2*ones(N,1),0) + diag(-1*ones(N-1,1),1));

counter = 1;
for j=h:h:1-h
  for i=h:h:1-h
  F(counter,1) =  f(i,j);
  counter = counter+1;
  endfor
endfor

counter2 = 1;
for j=0:h:1
  for i = 0:h:1
  usol(counter2,1) = u(i,j);
  counter2 = counter2+1;
  endfor
endfor

I = eye(N);

Ah2D = (kron(I,Ah) + kron(Ah,I));

% Computed solution
uh = Ah2D\F;

% Reshaping uh
U = reshape(uh,N,N);

% Now imposing boundary conditions
Utotal = zeros(N+2,N+2);
Utotal(:,1) = u(0:h:1,0);
Utotal(:,N+2) = u(0:h:1,1);
Utotal(1,:) = u(0,0:h:1);
Utotal(N+2,:) = u(1,0:h:1);
Utotal(2:N+1,2:N+1) = U;
Usol = reshape(usol,N+2,N+2);

% Plotting
[X,Y] = meshgrid(0:h:1, 0:h:1);
subplot(1,2,1)
surf(X,Y,Utotal)
subplot(1,2,2)
surf(X,Y,Usol)