% demo lin sys stable & unstable algs. Barnett 10/13/21
clear; rng(0);
A = ones(3,3) + 1e-14*rand(3,3)
c = [1;2;3];
b = A*c;
ct = inv(A)*b;                   % dangerous, may be unstable
norm(A*ct-b) / norm(b)           % rel resid terrible

ct = linsolve(A,b);              % use (backward) stable solver
norm(A*ct-b) / norm(b)           % rel resid great

norm(ct-c) / norm(c)             % rel err in soln vector: huge!

[U,S,V] = svd(A);
Y = diag(1./diag(S))*U';         % Y and V are factors of A
ct = V*(Y*b);                    % apply them
norm(A*ct-b) / norm(b)           % rel resid O(e_mach)
