u = N – u # Adjusting the initial guess and method to avoid overflow a_estimate_refined = fsolve(equation_to_solve, x0=0.001, args=(N, u))[0] # Calculate T = 1/a T = 1 / a_estimate_refined populations.append(T) return populations