I would like to fit a set of data (measured value vs time) with a model that contains a differential equation and an integral that has to be evaluated numericaly. The integral also contains Bessel functions (first and second kind) evaluated close to zero which I think does not help to solve the problem because of divergence.

The model is the following: $ $ \frac{dT\!\left(t\right)}{dt}=-\frac{16 c_0 Q}{\pi}\int_0^\infty \frac{\exp\left(-Qu^2t\right)}{u \left(J_0^2\left(u\rho\right)+Y_0^2\left(u\rho\right)\right)}T\!\left(t\right)^2 du $ $ $ u$ is just a dummy integration variable, and the parameters to be optimised are $ c_0, \rho$ and $ Q$ .

Here is the piece of code that I wrote, however without success:

(* Fit model *)  model = ParametricNDSolveValue[{T'[t] ==                                 NIntegrate[-((16 c Q)/Pi) Exp[-Q u^2 t]/                                            (BesselJ[0,u p]^2 + BesselY[0, u p]^2) T[t]^2,                                {u, 0, ∞}], T[0] = 1}, T, {t, 0, 2*^-6}, {c, Q, p}] (* fitting *)  fit = NonlinearModelFit[data,model[c, Q, p][t], {{c, 1*^15}, {Q, 1*^-6}, {p, 2*^-9}}, t] 

There are many errors, the first being that NIntegrate is evaluated at a non-numerical value.

Do you have suggestions to improve the code?

Here is a test data sample:

data = {{0, 1.05634}, {1/10000000, 0.752195}, {1/5000000, 0.495366},         {3/10000000, 0.484146}, {1/2500000, 0.342561}, {1/2000000, 0.43622},         {3/5000000, 0.235122}, {7/10000000, 0.128049}, {7.95*10^-7, 0.147317},         {9/10000000, 0.192805}, {1/1000000, 0.14561}, {1.1*10^-6, 0.32939},         {1.2*10^-6, 0.228659}, {1.3*10^-6, -0.0340244}, {1.4*10^-6, 0.0121951},         {1.5*10^-6, 0.131707}, {1.6*10^-6, 0.074878}, {1.7*10^-6, 0.30939},         {1.8*10^-6, 0.122561}, {1.9*10^-6, -0.00195122}, {2.005*10^-6, 0.0857317}}; 

**** EDIT ****

I have continued working on this problem. Now I first solve the integral, and I pass it to ParametricNDSolve. I do not think that evaluating Bessel functions close to 0 is a problem as the evaluation of the integral alone works well. The code is thus as follow:

f[c_?NumericQ, Q_?NumericQ, p_?NumericQ, t_?NumericQ] :=    NIntegrate[- ((16 c Q)/Pi) Exp[-Q u^2 t]/(BesselJ[0, u p]^2 + BesselY[0, u p]^2),   {u, 0, ∞}];  model = T /.ParametricNDSolve[{T'[t] == f[c, Q, p][t] T[t]^2, T[0] == 1},  T, {t, 0, 100}, {c, Q, p}];      fit = FindFit[data[[All, 2]],model[c, Q, p][t], {{c, 1*^15}, {Q, 1*^-6}, {p, 2*^-9}}, t]; 

Now I have the following error message: “ParametricNDSolve::pdvar Dependent variables {T,f[c,Q,p]} cannot depend on parameters {c,Q,p}”. I find it surprising because in principle one wants the dependent variables to depend on some parameters.