Big numbers problem while "odeint" numerical integration

Question

I'm having some computational problems with the following code:

import numpy as np
from numpy import arange
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from scipy.integrate import quad
import matplotlib as mpl
mpl.rcParams['agg.path.chunksize'] = 10000

# parameters

Ms = 100                                #GeV  Singlet Mass
Me = 0.511e-3                           #Gev Electron Mass
Mp = 1.22e19                            #GeV Planck Mass                                   
gs = 106.75                             # Entropy dof
H0 = 2.133*(0.7)*1e-42                  # GeV Hubble parameter (unused)
gx = 2                                  # WIMP's dof
g = 100                                 # total dof
sigmav=[1e-25,1e-11,1e-12]              # cross section's order of magnitude 

xi=1e-2
xe=1e2
npts=int(1e5)
x = np.linspace(xi, xe, npts)

def fMB(p,x,m):
    return np.exp(-x*np.sqrt(1+p*p/(m*m)))*p*p
def neq(x,m):
    return (gx/(2*np.pi*np.pi))*quad(fMB, 0, np.inf, args=(x,m))[0]
def neq_nr(x,m):
    return 2*(m**2/(2*np.pi*x))**(3/2)*np.exp(-x)
def stot(x):
    return (2*np.pi*np.pi/45)*gs*Ms*Ms*Ms/(x*x*x)
def Yeq(x,m):
    return neq(x,m)/stot(x)
Yeq2=np.vectorize(Yeq)
def Yeq_nr(x):
    return 0.145*(gx/gs)*(x)**(3/2)*np.exp(-x)
def Yeq_r(x):
    return 0.278*(3*gx/4)/gs
def Ytot(x):
    if np.any(x<=1):
        return Yeq_r(x)
    else:
        return Yeq_nr(x)

def eqd(yl,x,Ms,σv):
    '''
    Ms  [GeV]     : Singlet Mass
    σv: [1/GeV^2] : ⟨σv⟩
    '''                           
    H = 1.67*g**(1/2)*Ms**2/Mp
    dyl = -neq(x,Ms)*σv*(yl**2-Yeq(x,Ms)**2)/(x**(-2)*H*x*Yeq(x,Ms)) #occorre ancora dividere per Yeq_nr(x) oppure Yeq(x)
    return dyl


y0=1e-15
yl0 =  odeint( eqd, y0, x,args=(Ms,sigmav[0]), full_output=True)
yl1 =  odeint( eqd, y0, x,args=(Ms,sigmav[1]), full_output=True)
yl2 =  odeint( eqd, y0, x,args=(Ms,sigmav[2]), full_output=True)

fig = plt.figure(figsize=(11,8))
plt.loglog(x,yl0[0],  label = r'$\langle σ v\rangle = %s {\rm GeV}^{-2}$'%(sigmav[0]))
plt.loglog(x,yl1[0],  label = r'$\langle σ v\rangle = %s {\rm GeV}^{-2}$'%(sigmav[1]))
plt.loglog(x,yl2[0],  label = r'$\langle σ v\rangle = %s {\rm GeV}^{-2}$'%(sigmav[2]))
plt.loglog(x,Yeq_nr(x), '--', label = '$Y_{EQ}^{nr}$')
plt.loglog(x,Yeq2(x,Ms), '--', label = '$Y_{EQ}$')
plt.ylim(ymax=0.1,ymin=y0)
plt.xlim(xmax=xe,xmin=xi)
plt.xlabel('$x = m_χ/T$', size= 15)
plt.ylabel('$Y$', size= 15)
plt.title('$m_χ = %s$ GeV'%(Ms), size= 15)
plt.legend(loc='best',fontsize=12)
plt.grid(True)
plt.savefig('abundance.jpg',bbox_inches='tight', dpi=150)

In particular, as soon as I use little values of sigmav (ranging from 10^-12 to 10^-25) the solution is well displayed, but making use of bigger values (starting from 10^-11) I obtain problems and I guess is a order of magnitudes problem, but I don't know how to handle it!

Thanks to everyone!

Edit 1: I'm uploading a plot making use of three different values of sigmav and as you may see the bigger one (1e-10) is showing (I guess) precision problems plot_1