If X, Y, and Z are vectors of the x,y,z coordinates of your points, do this:
dX = X-mean(X);
dY = Y-mean(Y);
dZ = Z-mean(Z);
XX = sum(dX.^2); YY = sum((dY.^2; ZZ = sum(dZ.^2);
XY = sum(dX.*dY); XZ = sum(dX.*dZ); YZ = sum(dY.*dZ);
A = [XX,XY,XZ;
XY,YY,YZ;
XZ,YZ,ZZ];
[V,D] = eig(A); % <--Corrected
[~,k] = min(diag(D));
Then
V(1,k)*(x-mean(X))+V(2,k)*(y-mean(Y))+V(3,k)*(z-mean(Z)) = 0 % <--Corrected
where V(:,k) is the eigenvector of A with the smallest eigenvalue, and gives the equation for the best fit plane through your points. (Note: The lowercase x, y, and z here refer to coordinates of the equation of the plane.)
You can plot this plane by using an appropriate meshgrid of x,y points along with solving for z in the plane equation making use of the 'surf' command.
