clear

m=6;
A=HilbertMatrix(m); % hilb(m)
cond(A)
b=rand([m 1]);

x0=A\b;
x=InverseHilbertMatrix(m)*b;
norm(x-x0)/norm(x)

r=b-A*x0;
e=A\r;
x1=x0+e;
norm(x-x1)/norm(x)

r=b-A*x1;
e=A\r;
x2=x1+e;
norm(x-x2)/norm(x)