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);