
function phiu = Get_Phiu(ru,rl,K_approx,neta,shear_assumption)
    % Determines the critical value of Phiu given some constant ru, rl, neta
    % ru, rl = end fixity factors
    % K_approx = K for initial guess (use Newmark, 1949)
    % neta = EI/L^2 / kAG -- specify zero if no shear
    % shear_assumption = 0 for phi, 1 for dy/dx

    if isnan(rl) || isnan(ru) || isinf(rl) || isinf(ru)
        phiu = -163; % Invalid input
        return;
    end
    
    max_iter = 300;
    numIter = 0;
    converge = 0;
    
    % Initial guess
    phi_trial = (3.14159265358979 / K_approx);
    eps = phi_trial * 0.01;
    tol = 1e-12;
    retry = 0;
    
    while (numIter < max_iter && converge == 0)
        
        % Calculate omega as necessary
        if neta == 0  % No shear
            omega = 1;
        else % Yes shear
            if shear_assumption == 0 
                omega = sqrt(1 + neta * phi_trial ^ 2);
            else
                if (neta * phi_trial ^ 2) >= 1 
                    retry = 1;
                else
                    omega = 1 / sqrt(1 - neta * phi_trial ^ 2);
                end
            end
        end
        
        if retry == 1 
            break;
        end
        
        % Get slope at current point
        DEN_curr = beta_den(ru, rl, phi_trial, omega);
        if abs(DEN_curr) < tol 
            converge = 1;
            phiu = phi_trial;
        else
            DEN_step = beta_den(ru, rl, phi_trial + eps, omega);
            dDEN = DEN_step - DEN_curr;
            dDENdp = dDEN / eps;
            
            phi_next = phi_trial - DEN_curr / dDENdp;
            phi_trial = phi_next;
        end
    
        numIter = numIter + 1;
    end
    
    if (retry == 1 && K_approx < 2)  % Didn't converge and resulted in errors
        phiu = Get_Phiu(ru, rl, K_approx * 1.1, neta, shear_assumption);
    end
    
    if converge == 0 
        phiu = -100; % Didn't converge
    end
    if retry == 1 && K_approx >=2
        phiu = -103;
    end
end