function A=thirdsurf(x,y,Z)

[r,c]=size(x)
a = ones(r,1)
x1 = x
x2 = y
% transforming
for n = 1:r
    x3(n) = x(n)*x(n);
    x4(n) = x(n)*y(n);
    x5(n) = y(n)*y(n);
    x6(n) = x(n)*x(n)*x(n);
    x7(n) = x(n)*x(n)*y(n);
    x8(n) = x(n)*y(n)*y(n);
    x9(n) = y(n)*y(n)*y(n);
end

x3=x3'
x4=x4'
x5=x5'
x6=x6'
x7=x7'
x8=x8'
x9=x9'
X = [a,x1,x2,x3,x4,x5,x6,x7,x8,x9]
% computing the coefficients
A = inv(X'*X)*X'*Z

% R2 fit diagnosis
zg = X*A
SSD = (zg-Z)'*(zg-Z)
for n = 1:r
    zmean(n) = mean(Z)
end
zmean = zmean'
SSR = (zg-zmean)'* (zg-zmean)
pR = SSR/(SSR+SSD)

% F fit diagnosis
p = 9
f = (SSR/p)/(SSD/(r-p-1))
