Am I using poicalc correctly?

The solution returned using poicalc is too large by almost 1 order or magnitude. I expect the maximum z value of the surface to be 1, but instead I am getting 7.9976 (almost 8). I will post exactly what I did below.

I used the periodic, 2D function

And I set a = 1 to make things simpler

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); %sub in a = 1 into equation
u(x,y,a)

ans = 

%% Cartesian Space
% x = [0,.1,.2,...,X=1], y = [0,.1,.2,...,Y]
X = 1; Y = 1;
Lx = 99; Ly = 99;
dx = X/Lx; dy = Y/Ly;
x_arr = linspace(0,X,Lx+1);
y_arr = transpose(linspace(0,Y,Ly+1));
x_arr

x_arr = 1×100

0 0.0101 0.0202 0.0303 0.0404 0.0505 0.0606 0.0707 0.0808 0.0909 0.1010 0.1111 0.1212 0.1313 0.1414 0.1515 0.1616 0.1717 0.1818 0.1919 0.2020 0.2121 0.2222 0.2323 0.2424 0.2525 0.2626 0.2727 0.2828 0.2929

y_arr

y_arr = 100×1

0 0.0101 0.0202 0.0303 0.0404 0.0505 0.0606 0.0707 0.0808 0.0909

u_func(x,y,a) = u(x,y,a); % This is a term that we will be using to sub into u(x,y,a)

u(x,y,a) = subs(u(x,y,a),x,x_arr(1,:)); % sub in x 
u(x,y,a) = subs(u(x,y,a),y,y_arr(:,1)); % sub in y. I noticed that I could sub in a transposed array instead of using a meshgrid
u_xy = double(u(x,y,a)) % We need to convert this to a double for plotting purposes later on

u_xy = 100×100

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0016 0.0032 0.0047 0.0062 0.0077 0.0091 0.0105 0.0119 0.0132 0.0145 0.0158 0.0170 0.0182 0.0194 0.0206 0.0217 0.0228 0.0238 0.0248 0.0258 0.0267 0.0277 0.0285 0.0294 0.0302 0.0310 0.0317 0.0325 0.0331 0 0.0032 0.0063 0.0093 0.0123 0.0152 0.0180 0.0208 0.0235 0.0262 0.0288 0.0313 0.0337 0.0361 0.0385 0.0407 0.0429 0.0450 0.0471 0.0491 0.0511 0.0529 0.0547 0.0565 0.0582 0.0598 0.0613 0.0628 0.0642 0.0656 0 0.0047 0.0093 0.0138 0.0182 0.0225 0.0268 0.0309 0.0349 0.0389 0.0427 0.0464 0.0501 0.0536 0.0571 0.0604 0.0637 0.0669 0.0699 0.0729 0.0758 0.0786 0.0813 0.0839 0.0863 0.0887 0.0910 0.0933 0.0954 0.0974 0 0.0062 0.0123 0.0182 0.0241 0.0297 0.0353 0.0408 0.0461 0.0513 0.0563 0.0613 0.0661 0.0708 0.0753 0.0798 0.0841 0.0882 0.0923 0.0962 0.1000 0.1037 0.1072 0.1106 0.1139 0.1171 0.1201 0.1230 0.1258 0.1285 0 0.0077 0.0152 0.0225 0.0297 0.0368 0.0437 0.0504 0.0570 0.0634 0.0697 0.0758 0.0817 0.0875 0.0932 0.0986 0.1040 0.1091 0.1141 0.1190 0.1237 0.1282 0.1326 0.1368 0.1409 0.1448 0.1486 0.1522 0.1556 0.1589 0 0.0091 0.0180 0.0268 0.0353 0.0437 0.0519 0.0599 0.0677 0.0753 0.0827 0.0900 0.0970 0.1039 0.1106 0.1171 0.1234 0.1296 0.1355 0.1413 0.1468 0.1522 0.1574 0.1625 0.1673 0.1719 0.1764 0.1807 0.1848 0.1887 0 0.0105 0.0208 0.0309 0.0408 0.0504 0.0599 0.0691 0.0781 0.0869 0.0955 0.1038 0.1120 0.1199 0.1276 0.1352 0.1425 0.1495 0.1564 0.1630 0.1695 0.1757 0.1817 0.1875 0.1931 0.1984 0.2036 0.2085 0.2132 0.2178 0 0.0119 0.0235 0.0349 0.0461 0.0570 0.0677 0.0781 0.0883 0.0982 0.1079 0.1174 0.1266 0.1356 0.1443 0.1528 0.1610 0.1690 0.1768 0.1843 0.1916 0.1986 0.2054 0.2120 0.2183 0.2243 0.2301 0.2357 0.2411 0.2462 0 0.0132 0.0262 0.0389 0.0513 0.0634 0.0753 0.0869 0.0982 0.1093 0.1201 0.1306 0.1409 0.1508 0.1606 0.1700 0.1792 0.1881 0.1967 0.2051 0.2132 0.2210 0.2285 0.2358 0.2428 0.2496 0.2561 0.2623 0.2682 0.2739

%% Now we plot the solution
figure(1)
contour3(x_arr,y_arr,u_xy);
surface(x_arr,y_arr,u_xy);
title('u(x,y)');
xlabel('x = linspace(0,1,11)');ylabel('y = linspace(0,1,11)');

Then we set up the Poisson equation that we are trying to solve for in order to get back to the solution graphed above.

f_func = diff(u_func,x,2) + diff(u_func,y,2);
A = subs(f_func,'x',x_arr);
B = subs(A,'y', y_arr);
f = double(B);
size(f)

ans = 1×2

100 100

Now that we have f(x,y)< we will use poicalc on f(x,y) to get back to u(x,y).

%% Use poicalc function
% Now we run the poicalc function which takes in parameters (f,h1,h2,n1,n2)
h1 = dx; % gridspace_x
h2 = dy; % gridspace_y
n1 = 100; %number of points in x. supposed to be 100 but n1*n2 has to equal the number of rows in f
n2 = 1; % number of points in y. Supposed to be 100. 
           % n2 seems to dictate the number of local maxima and minima or
           % slopes
           % if I put n2 = 2, I get 2 slopes, and if n1 = n2 = 10, I get 10
           % slopes.
u_calc = poicalc(f,h1,h2,n1,n2); 
%% Now we plot it to examine how well it matches our first plot
figure(2)
contour3(y_arr,x_arr,abs(u_calc));
surface(y_arr,x_arr,abs(u_calc));
title('u(x,y) computed using poicalc function');
xlabel('x');ylabel('y');

max(max(abs(u_xy)))

ans = 0.9998

max(max(abs(u_calc)))

ans = 7.9976

Personally I don't have too much issue with the shape, my biggest issue is with the magnitude being off. The maximum value for u_xy in the first plot is very nearly 1, and the max for the solution obtained for u(x,y) here is very nearly 8. What exactly am I doing wrong with the poicalc function that is resulting in this disparity in magnitude?