Untitled

function [U c]=TriangularU(A, b)
[nr nc]=size(A);
W=[A b];

for k=1:nr-1
    for i=k+1:nr
        l(i,k)=W(i,k)/W(k,k);
        for j=k:nc+1
            W(i,j)=W(i,j)-l(i,k)*W(k,j);
        end
    end
end

U=W(:,1:nc);
c=W(:,nc+1);

end