math domain error in a python code involving fsolve and quad integrations

Im getting math domain error in this code. But the same function F does work in another program (I didn't use fsolve) without any domain error.

import numpy as np
from scipy.integrate import quad
from scipy.optimize import fsolve
from math import sqrt

def F(x, g, a, d, A, M):
  return A*(x**(1-g))/(1-g)+M*((x**(a+1))/(a+1)+(x**(d+1))/(d+1))

def C(q):
  return q+1

def intf1(x, rho, g, a, d, A, M):
  return (1/sqrt(F(rho, g, a, d, A, M)-F(x, g, a, d, A, M)))

def eq_to_solve(q, rho, g, a, d, A, M):
  return (quad(intf1,0,rho,args=(rho,g, a, d, A, M))[0]-quad(intf1,q,rho,args=(rho,g, a, d, A, M))[0]-(C(q)*q)/intf1(q, rho, g, a, d, A, M))

A = -4
M = 20
g = 0.6
a = 0.7
d = 2.1
rhomin = 10**(-6)
rhomax = 50
num_rho_values = 60
rholist = np.linspace(rhomin, rhomax, num=num_rho_values, endpoint=False)
lambdalist=[]
found_rholist = []
for rho in rholist:
    q_guess_list = np.linspace(1, rho, num = 10, endpoint=False)
    for q_guess in q_guess_list:
        q_val = fsolve(eq_to_solve, q_guess, args=(rho, g, a, d, A, M))
        if q_val > 0 and abs(eq_to_solve(q_val, rho, g, a, d, A, M)) < 10**(-6):
            lam = (C(q_val)*q)**2/(2*intf1(q_val, rho, g, a, d, A, M))
            print("found q= ", q_val)
            if lam > 0:
                print('found lam = ', lam)
                lambdalist.append(lam)
                found_rholist.append(rho)
            else:
                None
        else:
            None

I'm getting this error:

    62                                 q_guess_list = np.linspace(1, rho, num = 10, endpoint=False)
     63                                 for q_guess in q_guess_list:
---> 64                                         q_val = fsolve(eq_to_solve, q_guess, args=(rho, g, a, d, A, M))
     65                                         if q_val > 0 and abs(eq_to_solve(q_val, rho, g, a, d, A, M)) < 10**(-6):
     66                                                 lam = (C(q_val)*q)**2/(2*intf1(q_val, rho, g, a, d, A, M))

      6 
      7 def intf1(x, rho, g, a, d, A, M):
----> 8   return (1/sqrt(F(rho, g, a, d, A, M)-F(x, g, a, d, A, M)))
      9 
     10 def eq_to_solve(q, rho, g, a, d, A, M):

ValueError: math domain error