% Exploring accuracy in linear systems: A,b are the input data, x the output
% Barnett 4/11/13

A = [1 3;2 6+1e-14];
x = [2;3];
b = A*x 
xhat = A\b                % differs from known solution x badly

norm(x - xhat)            % forwards error bad, around 0.3
norm(x - xhat)/norm(x)    % fwd relative error also bad.

b = b + [1e-10;1e-11];    % perturb b a small bit
A\b                       % OMG, hugely wrong

cond(A)                   % this explains the sensitivity

norm(A*xhat - b)/norm(xhat)  % relative residual norm is close to eps_mach,
                             % as it must be for backwards-stable solve

A\randn(2,1)              % note typical RHS leads to huge ||x||
