Setup in python a function of Bessel functions that can be both integrated and differentiated

Question

I need to setup a user defined function ILTintegrand(u, t, r, b, a) for integrand = (j0(ur)*y0(ub) - y0(ur)*j0(ub)) / (j0(ub)*j0(ub) + y0(ub)*y0(ub)) * np.exp(-a*t*u**2) / u, where j0 and y0 are the Bessel functions of the first and second kind, respectively.

Then I need another user defined function ILTintegral(t, r, b, a) for the integral of ILTintegrand in u, bewteen 0 and infinity.

But then I also needed differentiating ILTintegral with respect to t and r, to check that a partial differential equation was satisfied (see this).

---

For the integration, I could do this using in ILTintegrand the Bessel functions from scipy with

from scipy.special import j0, y0

and then integrating in ILTintegral also with scipy with

from scipy import integrate
intg_pair = scipy.integrate.quad(ILTintegrand, 0, np.inf, args=(t, r, b, a))

But then I could not differentiate this wrt t or r.

---

When trying to integrate with sympy (so I could later use sympy.diff), using in ILTintegrand

from sympy import besselj, bessely

and then in ILTintegral with

intg = sympy.integrate(integrand, (u, 0, sym.oo))

the calculation was never completing.

---

What are possible ways of achieving my objective?

Related:

1. Differentiating a Bessel function with Lambda
2. Bessel function integral
3. How to calculate derivative and integral of the bessel functions in PYTHON

TL;DR

On my way of debugging the integration with sympy, I tried with shorter integrands, and changing the limits of integration to [0,1]. I found that

integrand = besselj(0, ur)
f = sym.integrate(integrand, (u, 0, 1))

works reasonably (though still slower than with scipy, I wouldn't know what happened if I had to evaluate this many times);

integrand = (besselj(0, ur) * bessely(0, ub) - bessely(0, ur) * besselj(0, ub))
f = sym.integrate(integrand, (u, 0, 1))

throws a long error report with two call stacks of many nested functions related to the computation of meijerg, like

...
ValueError: 0.302500000000000 is not an integer

During handling of the above exception, another exception occurred:

...
ValueError: expecting ints or fractions, got 0.302500000000000 and 1/2

and from

integrand = (besselj(0, ur) * bessely(0, ub) - bessely(0, ur) * besselj(0, ub)) \
    / ((besselj(0, ub))**2 + (bessely(0, ub))**2)
f = sym.integrate(integrand, (u, 0, 1))

towards the final destination, none of the integrations finished.