digits(100); rw = []; vw = [];
for n = 3:20,
  A = hilb(n); b = ones(n,1);
  x = A\b; % Solve LSE with std precision
  rw = [rw;n,norm(A*x-b)/norm(x)];
  xe = vpa(A)\vpa(b); % high precision solve
  vw = [vw;n,norm(double(x-xe))/norm(double(xe))];
end
figure('name','Gaussian elimination for Hilbert matrix');
plot(rw(:,1),rw(:,2),'b-',vw(:,1),vw(:,2),'r--');
xlabel('n','fontsize',14);

