%
% Solve the linear system Ux = b using back substitution.
%
% Input:  Upper triangular matrix U
%         vector b
%
% Output: vector x satisfying Ux = b
%
function x = backsub(U, b)
    
    % Determine the number of rows/columns in U
    [n, m] = size(U);
    
    % Ensure the matrix U is square
    if n ~= m
        error('Square matrix expected.');
    end
    
    % Create a column vector x of the appropriate size
    x = zeros(n, 1);
    
    % Solve for x, starting with x(n) and working toward x(1)
    for i = n : -1 : 1
        sum = 0;
        for j = i + 1 : n
            sum = sum + U(i, j) * x(j);
        end
        x(i) = (b(i) - sum)/U(i, i);
    end
end