Question: Originally, I had the following equation: 2mgv×sin(\alpha) = CdA×\rho(v^2 + v_{wind}^2 + 2vv_{wind}cos(\phi))^(3/2)  which I could express as the following non-linear equation: (K × v)^(2/3) = v^2 + v_{wind} + 2vv_{wind}cos(\phi)  To solve this, I need to use a numerical approach. I tried writing a code for this in Python using fsolve from scipy.optimize. However, the results I got are not too promising. What else should I try? Should I use a different approach/package or just my code needs to be improved? I also experienced that the result is highly dependent on v_initial_guess. Please note, that I would consider myself as a beginner in programming. I also tried writing a code to solve the equation for v using the Newton-Raphson method but I wasnt too successful. Here is my code: import numpy as np from scipy.optimize import fsolve m = 80 g = 9.81 alpha = np.radians(10) #incline CdA = 0.65 rho = 1.225 v_w = 10 phi = np.radians(30) #wind angle with the direction of motion sin_alpha = np.sin(alpha) cos_phi = np.cos(phi) def equation(v): K = ((m * g * sin_alpha) / ((CdA * rho))**(2/3)) return K * v**(2/3) - v**2 - 2*v*v_w*cos_phi - v_w**2 v_initial_guess = 30 v_solution = fsolve(equation, v_initial_guess, xtol=1e-3) print("v:", v_solution[0])type here  EDIT: This is what my code looks like now:  import numpy as np from scipy.optimize import fsolve import matplotlib.pyplot as plt m = 80 g = 9.81 alpha = np.radians(2) # incline CdA = 0.321 rho = 1.22 v_w = 5 phi = np.radians(180) # wind angle with the direction of motion sin_alpha = np.sin(alpha) cos_phi = np.cos(phi) def lhs(v): return m * g * v * sin_alpha def rhs(v): return 0.5 * CdA * rho * (v**2 + v_w**2 + 2*v*v_w*cos_phi)**(3) def difference(v): return lhs(v) - rhs(v) # fsolve to find the intersection v_initial_guess = 8 v_intersection = fsolve(difference, v_initial_guess)[0] v_values = np.linspace(0.1, 50, 500) lhs_values = lhs(v_values) rhs_values = rhs(v_values) plt.figure(figsize=(10, 6)) plt.plot(v_values, lhs_values, label='$2mgv\\sin(\\alpha)$', color='blue') plt.plot(v_values, rhs_values, label='$CdA\\rho(v^2 + v_{wind}^2 + 2vv_{wind}\\cos(\\phi))^{3/2}$', color='red') plt.xlabel('Velocity (v)') plt.xlim(0, 20) plt.title('LHS and RHS vs. Velocity') plt.legend() plt.grid(True) plt.ylim(0, 2000) plt.show() print(f"The intersection occurs at v = {v_intersection:.2f} m/s") P_grav_check = m *g * sin_alpha * v_intersection P_air_check = 0.5 * CdA * rho * (v_intersection ** 2 + v_w ** 2 + 2 * v_intersection * v_w * cos_phi) ** (3) print(P_grav_check) print(P_air_check) 

Answer: The solution, using the secant method is the following: import numpy as np import matplotlib.pyplot as plt m = 60 g = 9.81 alpha = np.radians(2) # incline CdA = 0.321 rho = 1.225 v_wind = 4 phi = np.radians(60) # wind angle with the direction of motion mu = 0.005 cos_phi = np.cos(phi) sin_alpha = np.sin(alpha) lhs = m * g * sin_alpha # gravitational force def rhs(v): return 0.5 * CdA * rho * (v + v_wind *cos_phi)**2 # drag force def difference(v): return lhs - rhs(v) def secant_method(func, v0, v1, tolerance=1e-6, max_iterations=1000): for _ in range(max_iterations): f_v0 = func(v0) f_v1 = func(v1) v_new = v1 - f_v1 * (v1 - v0) / (f_v1 - f_v0) if abs(v_new - v1) < tolerance: return v_new v0, v1 = v1, v_new raise ValueError("Secant method did not converge") v_initial_guess_1 = 15 v_initial_guess_2 = 25 v_intersection = secant_method(difference, v_initial_guess_1, v_initial_guess_2) v_values = np.linspace(0.1, 50, 500) rhs_values = rhs(v_values) P_grav_check = m * g * sin_alpha * v_intersection P_air_check = 0.5 * CdA * rho * (v_intersection + v_wind * cos_phi)**2 * v_intersection P_friction = mu * m * g * np.cos(alpha) * v_intersection alpha_degree = np.degrees(alpha) alpha_percentage = np.tan(alpha) *100 print(f"\nopposing component of v_wind: ",v_wind *cos_phi, "m/s") print(f"The intersection occurs at v = {v_intersection:.4f} m/s") print(f" at a {alpha_degree} degrees incline which equals to {alpha_percentage:.2f} %") print(f"\nP_grav: {P_grav_check:.2f} W") print(f"P_air_check: {P_air_check:.2f} W") print(f"P_friction: {P_friction:.2f} W") print(f"P = {P_grav_check + P_air_check + P_friction}") plt.figure(figsize=(10, 6)) plt.plot(v_values, [lhs]*len(v_values), label='$mg\\sin(\\alpha)$', color='blue') plt.plot(v_values, rhs_values, label='$0.5CdA\\rho ({v + v_{wind}\\cos(\\phi)})^2$', color='red') plt.xlabel('Velocity (v)') plt.xlim(0, 20) plt.title('Gravitational vs. Drag Force') plt.legend() plt.grid(True) plt.ylim(0, 200) plt.show()