function x = backsub(A,b)
%
%  Solve the upper-triangular system A*x = b
%  by backward substitution
%
   n = size(A,1);  % the row dimension of A, assuming A is square
   x = zeros(n,1); % create a column vector of dimension n
   x(n) = b(n)/A(n,n);   
   for k = (n-1):-1:1
      x(k) = (b(k) - A(k,k+1:n)*x(k+1:n))/A(k,k);
   end