function A=secondsurf(x,y,Z)


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

x3=x3'
x4=x4'
x5=x5'
X = [a,x1,x2,x3,x4,x5]
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 = 5
f = (SSR/p)/(SSD/(r-p-1))

