function Pu = Iterate_Pu(E0,Fy,nu,L,I,A,ru0,rl0,with_shear,kappa,shear_assumption,bool_buckle)
    % Determines the value Pu with considering the tangent modulus reduction (Yura & Helwig, 2001) for a semi-rigid column
    % shear_assumption = 0 for phi, 1 for dy/dx
    
    max_iters = 100;
    num_iters = 0;
    converge = 0;
    tol = 0.001;
    E = E0;
    ru = ru0;
    rl = rl0;
    Py = Fy * A;
    Pu = 0;
    e_val = 0;
        
    % Get the beam rotational restraints
    RBu =0;
    RBl =0;
    if ru0 ~= 1 
        RBu = 3 * E0 * I / L * ru / (1 - ru);
    end
    if rl0 ~= 1 
        RBl = 3 * E0 * I / L * rl / (1 - rl);
    end
    
    while (num_iters < max_iters && converge == 0)
        
        % Update end fixity factor with new E
        if ru0 == 1 
            ru = 1;
        elseif ru0 == 0 
            ru = 0;
        else
            ru = 1 / (1 + 3 * E * I / (L * RBu));
        end
        if rl0 == 1 
            rl = 1;
        elseif rl0 == 0 
            rl = 0;
        else
            rl = 1 / (1 + 3 * E * I / (L * RBl));
        end
        
        % Update G
        if bool_buckle == 1 
            e_val = 0.5 * (E0 / E - 1);
        end
        G = E0 / (2 * (1 + nu) + 3 * e_val);
        kAG = kappa * A * G;
        
        % Solve phiu function with the updated end fixity factors
        K_approx = Get_K(ru, rl);
        if with_shear == 1 
            neta = E * I / L ^ 2 / (kAG);
        else
            neta = 0;
        end
        phiu = Get_Phiu(ru, rl, K_approx, neta, shear_assumption);
        
        % Calculate the ultimate load load, Pu, for the given E
        Fe0 = phiu ^ 2 * (E0 * I / L ^ 2) / A;
        new_Pu = 0.85 * Py * ((10 ^ (-1 / 7.38)) ^ (Fy / Fe0));
        if abs(new_Pu - Pu) < tol 
            converge = 1;
        end
        Pu = new_Pu;
        
        % Update E
        if (Pu / Py < 0.3333) 
            E = E0;
        else
            E = E0 * (-7.38 * Pu / Py * log(Pu / Py / 0.85) / log(10)); % VBA uses lnx
        end
        
        num_iters = num_iters + 1;
    end
    
    if converge == 0 
        Pu = -100; % Didn't converge error
    end
end
