You can try using matlab's 'svd' function. It will find a set of coefficients A11, A22, A33, ...,A0 for which the sum of the squares of the left sides of your quadratic expression for each of your (x,y,z) points is a minimum subject to the constraint that the sum of the squares of the coefficients be one. There are many ways of defining a "best" fit but this is one definition that has an easy solution.
Let X, Y, Z be column vectors for your given points. That is, the n-th elements of the three vectors are, respectively, the x, y, and z coordinates for the n-th point. Then form the matrix M and use 'svd' to do a singular value decomposition of M:
M = [X.^2,Y.^2,Z.^2,X.*Y,X.*Z,Y.*Z,X,Y,Z,ones(size(X))];
[~,S,V] = svd(M,0);
A = V(:,10);
The vector A which is the rightmost column of V will contain the desired list of coefficients, A11 A22, etc. It is a normalized eigenvector and the sum of its squares will accordingly be one. The smallest singular value that occurs at the lower right side of the main diagonal of the 10-by-10 matrix S, S(10,10), can be considered a measure of the quality of fit.
Best Answer