function x = bs1(U,c)
% BS1  backsubstitution with some improvements:
%        (i)   replace inner loop with colon notation
%        (ii)  stop if divide-by-zero
%        (iii) check size compatibility of inputs
%        (iv)  check if U is upper triangular
% usage:   >> x = bs1(U,c)     % should be equivalent to  x = U\c

n = size(U,1);
if size(U,2) ~= n, error('matrix U not square'), end
if size(c,2) ~= 1, error('vector c is not a column vector'), end
if size(c,1) ~= n, error('vector c not compatible size with U'), end

if max(max(abs(tril(U,-1)))) ~= 0, error('matrix U not upper triangular'), end

x = zeros(n,1);

scale = max(max(abs(U)));       % size of largest entry in U
if abs(U(n,n))/scale < 1e-14
  error('U(n,n) is zero; divide-by-zero!')
end

x(n) = c(n) / U(n,n);

for j = n-1:-1:1                % j is the row number
  sum = U(j,j+1:n) * x(j+1:n);  % replaces entire inner loop in bs0.m
  if abs(U(j,j))/scale < 1e-14
    error(sprintf('U(%d,%d) is zero; divide-by-zero!',j,j))
  end
  x(j) = (c(j) - sum) / U(j,j); % U(j,j) is diagonal element in row j
end

