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'*d)*Z

contour(X1(1,:),X2(:,1),S',10); set(gca(),"auto_clear","off");
contour(X1(1,:),X2(:,1),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:");
