ifft2 not returning expected plots

To start with, I am trying to solve the Poisson Equation using fft2 and ifft2. I am starting with the solution, with a=1, taking the Laplacian of u(x,y) to get and then trying to use the fft2 and ifft2 functions to work back to u(x,y). A more in-depth description of the problem and process can be found in the top answer found here. The problem comes at the 4th step in the linked answer when I take ifft2. The shape of the plot that I obtain for u(x,y) is wrong and the values are way too large. I will post my code below so that anyone can reproduce the problem as well as the plots of what I have obtained in the steps leading up to the problem.

syms u(x,y,a)
u(x,y,a) = 2^(4*a) * x^a * (1-x)^a * y^a * (1-y)^a;
u(x,y,a) = subs(u(x,y,a),a,1);
X = 1;Y=1;
Lx = 10; Ly = 10;
dx = X/Lx; dy = Y/Ly;
x_arr = meshgrid(0:dx:X,0:dx:X);
y_arr = transpose(meshgrid(0:dy:Y,0:dy:Y));

u_func(x,y,a) = u(x,y,a);

u(x,y,a) = subs(u(x,y,a),x,x_arr(1,:));
u(x,y,a) = subs(u(x,y,a),y,y_arr(:,1));

u_xy = double(u(x,y,a));

figure(1)
contour3(x_arr,y_arr,u(x,y,a));
surface(x_arr,y_arr,u(x,y,a));
title('u(x,y)');
xlabel('x = linspace(0,1,11)');ylabel('y = linspace(0,1,11)');

f_func = diff(u_func,x,2) + diff(u_func,y,2);
f_func = subs(f_func,x,x_arr(1,:));
f_func = subs(f_func,y,y_arr(:,1));
 
M = 10;
N = 10;
%m = meshgrid(linspace(0,M,M+1),linspace(0,N,N+1));
%n = transpose(meshgrid(linspace(0,N,N+1),linspace(0,N,N+1)));
 
P = M; Q = N;
Lp = P; Lq = Q;
dp = Lp/P;dq = Lq/Q;
p = meshgrid(0:dp:P,0:dp:P);
q = transpose(meshgrid(0:dq:Q,0:dq:Q));
 
f = double(f_func);
fhat = fft2(f,P+1,Q+1); %m = n = 10
 
km_axis = (-1j * 2 * pi/length(x_arr)).*p;
ln_axis = transpose((-1j * 2 * pi/length(y_arr)).*q);

figure(2)
contour3(real(km_axis) + imag(km_axis),real(ln_axis) + imag(ln_axis),abs(fftshift(fhat)));
surface(real(km_axis) + imag(km_axis),real(ln_axis) + imag(ln_axis),abs(fftshift(fhat)));
title('f(x,y) = $\nabla^2 u(x,y)$','Interpreter','LaTex')
xlabel('$k_{m}$ = $-j2 \pi p/Lx$','Interpreter','LaTex');ylabel('$l_{n}$ = $-j2 \pi q/Ly$','Interpreter','Latex');

uhat = -fhat./ (km_axis.^2 + ln_axis.^2);
mask = isinf(uhat);
uhat(mask) = 0;%assume that the values are 0 along the boundary
u = ifft2(uhat);

figure(3)
contour3(x_arr,y_arr,abs(u));
surface(x_arr,y_arr,abs(u));
title('$u(x,y) = ifft2\frac{\hat{f}}{(k_m^2 + l_n^2)}$','Interpreter','LaTex')
xlabel('x');ylabel('y');

I'm expecting my plot for figure 3 to match my plot for figure 1 in both shape and values. What I keep getting instead is a 3D parabola with values which are way too large and I am not certain as to why this is the case. I have been hung up on this for at least a month now and don't understand what I could be doing wrong whether it be in my methodology or if I am making a mistake somewhere in the code.

One thing I have noted is that if I change p and q to be centered at 0 such that:

p = meshgrid(-P/2:dp:P/2,-P/2:dp:P/2);
q = transpose(meshgrid(-Q/2:dq:Q/2,-Q/2:dq:Q/2));

I get a result for figure 3 that is much closer to figure 1 though I do not know why this is the case. Regardless, I cannot seem to figure out what is wrong with my approach whether it be in the methodology or in my usage of the functions.