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

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

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

contour(X1,X2,S); hold on;
contour(X1,X2,h,[0 0]); xlabel('x1'); ylabel('x2');
plot([0;d(1)],[0;d(2)],'r',[0;dp(1)],[0;dp(2)],'g',...
     [d(1),dp(1)],[d(2),dp(2)],'k:');
