function x = backsolve(A,b)
% solve the upper triangular linear system Ax = b
n = length(b);
if any(size(A) ~= [n n])
    error('lengths not compatible')
end
x = zeros(n,1); % column vector, otherwise get error
for j=n:-1:1
    rhs = b(j) - A(j,j+1:n)*x(j+1:n); % dot product of jth row of A 
                                      % starting to the right of the
                                      % diagonal with the part of x that is
                                      % already known
    x(j) = rhs/A(j,j);
end