x=0:0.1:10;
[X2,X3]=meshgrid(x,x);
X1=X2;
S=(X1-5).^2+(X2-5).^2+(X3-3)^2;

M=[1; -1; 0]
[Q,R]=qr(M)
Z=Q(:,2:3)
H=[2 0 0;0 2 0;0 0 2]
Hp=Z'*H*Z

d=2*[2;3;1]
dp=Z*(Z'*d)

mesh(X1,X2,X3); hold on;
xlabel('x1'); ylabel('x2'); zlabel('x3');
plot3([0;d(1)],[0;d(2)],[0;d(3)],'r',[0;dp(1)],[0;dp(2)],[0;dp(3)],'g',...
     [d(1),dp(1)],[d(2),dp(2)],[d(3),dp(3)],'k:');
