function x = bksolve(A, b)
% Solves Ax = b for x, assuming that A is a nonsingular upper 
% triangular matrix.

n = length(b);
x = b;

% add tests for A(i,i) close to 0!
x(n) = b(n) / A(n,n);
for i = n-1:-1:1
    x(i) = (b(i) - A(i,i+1:n)*x(i+1:n)) / A(i,i);
end

return

% Try
%    A = [2,2,3,4;0,5,6,7;0,0,8,9;0,0,0,10];
%    b = [20,34,25,10]';
%    bksolve(A,b)
%
% should give
%  ans =
%
%       2
%       3
%       2
%       1

