2D FFT square function

I am trying to perform the 2D FFT on the square function. I have defined the square function using:

if true
deltax = L/N;
deltay = L/N;
x = [-L/2:deltax:L/2-deltax];
y = [-L/2:deltay:L/2-deltay];
[X,Y]=meshgrid(x,y);
func = double(abs(X)<= b/2).*double(abs(Y)<= b/2);
figure;
surf(X,Y,func);
end

This does give the expected answer, namely:

However, when I go further and try and do the FFT using the following code:

if true
  % code
end
deltafx = 1/L;
deltafy = 1/L;
fx = [-1/(2*deltax):deltafx:1/(2*deltax)-deltafx];
fy = [-1/(2*deltay):deltafy:1/(2*deltay)-deltafy];
[FX,FY]=meshgrid(fx,fy);
funcfourier=deltax*deltay*fft2(func);
funcfourier=fftshift(funcfourier);
absfuncfourier=abs(funcfourier);
surf(FX,FY,absfuncfourier);

gives:

which is not what I expected as I cannot read any amplitudes of that plot. I do not really know what went wrong here.

Any suggestion of how to show this in 3D is appreciated. Thanks

[jyoti mundra]

Please let me know what is L and N in your code?