function p = ModelOptimize(x,y,z)
x = x(:);
y = y(:);
z = z(:);
l = ones(numel(x),1);
o = zeros(numel(x),1);
% data
X = [ x.^2, y.^2, z.^2, x.*y, x.*z, y.*z, x, y, z, l ];
% gradients
dx = [2*x, o, o, y, z, o, l, o, o, o];
dy = [o, 2*y, o, x, o, z, o, l, o, o];
dz = [o, o, 2*z, o, x, y, o, o, l, o];
dX = [dx; dy; dz];
% scatter matrices
A = X'*X;
B = dX'*dX;
[V,D] = eig(A,B);
[~,ix] = min(abs(diag(D)));
p = V(:,ix);
