3D quadratic fitting to 27 voxel points, extend Matlab's findpeak function to 3rd dimention

the function below from Matlab's findpeak, fits a 2D quadtratic function u(x,y) to 9 pixels, I would like to extend that to 3D. i.e. u(x,y,z) for a 3x3x3 matrix.

-----------------------------------------------------------------

% simplified findpeak function below

[x_max,y_max] = simple_findpeak(u)

% fit a 2nd order polynomial to 9 points

% u is a 3x3 matrix

u = u(:);

x = [-1 -1 -1 0 0 0 1 1 1]';

y = [-1 0 1 -1 0 1 -1 0 1]';

% u(x,y) = A(1) + A(2)*x + A(3)*y + A(4)*x*y + A(5)*x^2 + A(6)*y^2

X = [ones(9,1), x, y, x.*y, x.^2, y.^2];

% u = X*A

A = X\u;

% get absolute maximum, where du/dx = du/dy = 0

x_max = (-A(3)*A(4)+2*A(6)*A(2)) / (A(4)^2-4*A(5)*A(6));

y_max = -1 / ( A(4)^2-4*A(5)*A(6))*(A(4)*A(2)-2*A(5)*A(3));