ValueError: Residuals are not finite in the initial point. (Python)

Question

I need to fit a graph using a complicated function including numerical integration. If I do this using a simpler function (which is also an integral), this method works, but for the function I'm trying, it's not working. Although the program is running, it is making a very bad fit. You can see that the adjustment the software makes is a constant line y = 0 ... this is not expected.

The code is:

import numpy as np
from scipy.integrate import quad
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import pandas as pd

dados=pd.read_excel("C:/Users/Samaung/Documents/Henrique/DadosExcel/dadosbalanca.xlsx")
t=dados["t"].values
m=dados["m"].values

mt=3.75*10**-5 # kg massa total empregada da amostra
rho=5961 #kg/m**3 densidade da amostra
v0=mt/rho #m**3 volume da amostra
uim=0.453 #A*m**2 momento magnético ímã
g=9.82 #m/s**2 aceleração da gravidade na Terra
z0=0.0063 #m distância centro ímã parte inferior da amostra
z1=0.0093 #m distância centro ímã partesuperior da amostra
R=0.0049 #m raio do porta amostra cilíndrico
dis=1/z0**4-1/z1**4-(z0**2+R**2)/(3*(z0**2+R**2)**3)+(z1**2+R**2)/(3*(z0**2+R**2)**3)
Msv=470644.784 #A/m magnetização volumétrica
Dm=17.6*10**-9 #m diâmetro médio
sigmad=0.13 #fator de dispersão adimensional
u0=1.256*10**-6 #N/A**2 permeabilidade magnética
H=70823.95 #A/m campo aplicado
Kb=1.3807*10**-23 #mm*2*kg/(s*K) constante de Boltzman
T=300 #K temperatura
Hk_=2/(u0*Msv) #campo de anisotropia dividido por Kef
t0=10**-9 #s fator de frequência
h_=H/Hk_
v1=6*Kb*T

def f(t,Kef):
    Multiplicador=3*u0*v0*uim/(32*np.pi*g)*dis*Msv*np.exp(2*sigmad**2)/(Dm*np.sqrt(2*np.pi)*sigmad)
    Integrando1=lambda D: np.exp(-(np.log(D/Dm))**2/(2*sigmad**2))
    Integrando2= lambda D: 1/np.tanh(u0*np.pi*Msv*H*D**3/v1)-v1/(u0*np.pi*Msv*H*D**3)
    Integrando3=lambda D: 1-np.exp(-t*((1-(h_/Kef)**2)*((1-h_/Kef)*np.exp(-Kef*np.pi*D**3/v1*(1-h_/Kef)**2)+(1+h_/Kef)*np.exp(-Kef*np.pi*D**3/v1*(1+h_/Kef)**2)))/(t0))
    fn=lambda D: Integrando1(D)*Integrando2(D)*Integrando3(D)
    ma=np.asarray(quad(fn,Dm/10,10*Dm)[0])
    return Multiplicador*ma
func=np.vectorize(f)

parametro, erro=curve_fit(func,t,m, p0=(139178.61676287447,), bounds=(0.1, np.inf))
print("Kef =", parametro[0])
xFit=np.arange(0.1,27105.06,1)
y = func(xFit, *parametro)
plt.plot(t, m, 'o')
plt.plot(xFit, y, '-')
plt.show()

My data are: dados= https://github.com/Henrique0501/Dadosbalanca/blob/main/dadosbalanca.xlsx

Could you help me to fix this?