% Wilkinson's counterexample: instability of Gaussian elimination
% despite column pivoting
res = [];
for n=10:10:200
  % Build Wilkinson matrix
  A = [tril(-ones(n,n-1))+2*[eye(n-1);zeros(1,n-1)],ones(n,1)];
  % imposed solution
  x = ((-1).^(1:n))'; b = A*x;
  % Numerical solution
  xh = A\b;
  % Residual
  r = b-A*xh;
  % Relative forward error
  relerr = norm(xh-x)/norm(x);
  % Relative norm of residual
  relresnorm = norm(r)/norm(xh);
  res = [res; n relerr relresnorm];
end
semilogy(res(:,1),res(:,2),'m-*',res(:,1),res(:,3),'r-+');
set(gca,'fontsize',14);
title('Numerical error: LSE with Wilkinson matrix');
xlabel('matrix size n','fontsize',14);
ylabel('relative errors/norms','fontsize',14);
legend('relative forward error','relative residual norm',...
       'location','best');
