Python (using SciPy) is unable to solve for a variable in an equation

I am having a problem solving for the variable x_Prime using the scipy.optimize.fsolve() function. Below is what I have:

import math
import scipy.optimize as so

c_val = 100
f_NFW_val = math.log(1+c_val)-c_val/(1+c_val)

psi_prime = 3.0608421399604424

eqn = lambda x_Prime: psi_prime  -  1/f_NFW_val * ( math.log(1+c_val*x_Prime)/x_Prime - c_val/(1+c_val*x_Prime) )
sol = so.fsolve(eqn, 1)[0]

I get the error:

It seems like there is no solution for psi_prime = 3.0608421399604424. However, if I try plotting the curve using the following code:

c_val = 100
f_NFW_val = math.log(1+c_val)-c_val/(1+c_val)

psi_prime = 3.0608421399604424

xs = np.linspace(0.1633804055843348, 1, 100)
plt.plot(xs, 1/f_NFW_val * ( np.log(1+c_val*xs)/xs - c_val/(1+c_val*xs) ))
plt.plot(xs, psi_prime*np.ones(np.size(xs)), color="grey", linestyle="dashed")
plt.show()

So, there should be a solution. Can someone tell me why fsolve() is not able to solve the equation?

The domain error is caused by an invalid argument of the logarithm. You should use numpy.log to handle cases like this (it will return NaN instead of raising an exception).

The second point is that you use an approach that is likely to diverge and stuck somewhere around the initial point. You should try different methods and analyze the function graph to come up with a way to find the root numerically, there is just no universal method to do so. In your case a default method from scipy.root_scalar should work just fine.

eqn = lambda x_Prime: psi_prime  -  1/f_NFW_val * ( np.log(1+c_val*x_Prime)/x_Prime - c_val/(1+c_val*x_Prime) )
sol = so.root_scalar(eqn, bracket=(0, 2))
print(sol)

It converges at 0.17997695477705547

root_scalar documentation