A=[1 -1 1;
1 -0.5 .25;
1 0 0;
1 .5 .25;
1 1 1];
b=[1 .5 0 .5 2.0]';
v=[sqrt(5)+1 1 1 1 1]';
v1=A(:,1)+norm(A(:,1),2)*[1 0 0 0 0]';
H1 = eye(5)- 2 *v*v'/(v'*v);
A2=H1*A
% Another way without computing H1
A22=A - 2*(v*(v'*A)/(v'*v));
v2=[0, A2(2:5,2)']'-norm(A2(2:5,2),2)*[0 1 0 0 0]';
H2= eye(5)- 2 *v2*v2'/(v2'*v2);
A3=H2*A2;
v3=[0,0,A3(3:5,3)']'+ norm(A3(3:5,3),2)*[0 0 1 0 0]';
H3= eye(5)- 2 *v3*v3'/(v3'*v3);
A4=H3*A3;
B=H3*H2*H1*b;
x=[A4(1:3,:)] B(1:3)