function [x] = backsub(u,b);
% Usage: [x] = backsub(u,b);
% Solve upper triangular system ux = b by back substitution
% matrix is accessed column-wise, and information from previously
% calculated elements of x is integrated progressively by modifying b.
% Input:
% u = upper triangular matrix
% b = right hand side
% Output:
% x = solution vector

[n,m] = size(u); % Find size of system

for j = n:-1:1 % Loop backwards over columns
    if u(j,j)==0 % Stop if matrix is singular
       disp('Matrix is singular'); 
       break
    end
    x(j) = b(j)/u(j,j);
    b(1:j-1) = b(1:j-1) - u(1:j-1,j)*x(j); % Update right hand side
end  
x=x'; % Transform x from row to column vector