function x=BackSubstitutionSAXPY(U,b)
% BACKSUBSTITUTIONSAXPY solves linear system by backsubstitution 
%   x=BackSubstitutionSAXPY(U,b) solves Ux=b by backsubstitution by
%   modifying the right hand side (SAXPY variant)

n=length(b);
for i=n:-1:1
  x(i)=b(i)/U(i,i);
  b(1:i-1)=b(1:i-1)-x(i)*U(1:i-1,i);
end 
x=x(:);