function X = uptrbk(A, B)
    % initialization
    [N N] = size(A);
    X = zeros(N, 1);
    C = zeros(1, N+1);
    
    % construct augmented matrix
    Aug = [A B];
    
    disp('Aug before elimination');
    disp(Aug);
    
    for p=1:N-1
        % partial pivoting
        [Y,j] = max(abs(Aug(p:N, p)));
        C = Aug(p,:);
        Aug(p,:) = Aug(j+p-1,:);
        Aug(j+p-1,:) = C;
        
        disp('Aug after pivoting');
        disp(Aug);
        
        if Aug(p, p) == 0
            'A was singular';
            break
        end
        
        % Elimination
        for k = p+1:N
            m = Aug(k, p) / Aug(p,p);
            Aug(k, p:N+1) = Aug(k, p:N+1) - m*Aug(p,p:N+1);
        end
        
        disp(['Aug after ', num2str(p), ' elimination']);
        disp(Aug);
    end
    
    X = backsub(Aug(1:N, 1:N), Aug(1:N, N+1));
end

function X = backsub(A, B)
    n = length(B);
    X = zeros(n, 1);
    X(n) = B(n)/A(n, n);
    
    for k = n-1:-1:1
        X(k) = (B(k) - A(k, k+1:n) * X(k+1:n)) / A(k,k);
    end
end