I often encounter the integrals in the following form:

$\int_0^\infty{\rm Bessel}(ax)\cdot{\rm Bessel}(bx)\cdot f(cx)dx$,

where Bessel can be $J$, $N$, $H^{(1)}$, $H^{(2)}$, $I$, or $K$; and $f(x)$ can be $\sin(x)$, $e^x$, etc. For example,

$\int_0^\infty K_\nu(ax)I_\nu(bx)\cos(cx)dx=\frac{1}{2\sqrt{ab}}Q_{\nu-1/2}(\frac{a^2+b^2+c^2}{2ab})$,$\qquad{\rm Re}(a)>|{\rm Re}(b)|$, $c>0$, ${\rm Re}(\nu)>-1/2$

I already know the result of this integral because it is in Gradshtein & Ryzhik's book. Sometimes the integration is with respect to the order of the Bessel function.

Both Mathematica 8 and Maple 15 cannot do this kind of integrals. When the integral involves two Bessel functions or two other special functions, Mathematica and Maple usually cannot do even if the integral has a closed-form result. My questions are as follows:

1. Is there any general theory about how to do this kind of integrals? (I wonder how the authors of that book did the above integral.)

2. I know there is a Mathematica package "Holonomic Functions" (http://www.risc.jku.at/research/combinat/software/HolonomicFunctions/). It seems that this package can help, but it does not seem very straightforward to obtain the final result. This package can verify the integrals with already known results, but can it do new integrals as above? Are there any better ways to deal with these integrals by computer?

I often encounter the integrals in the following form:

$\int_0^\infty{\rm Bessel}(ax)\cdot{\rm Bessel}(bx)\cdot f(cx)dx$,

where Bessel can be $J$, $N$, $H^{(1)}$, $H^{(2)}$, $I$, or $K$; and $f(x)$ can be $\sin(x)$, $e^x$, etc. For example,

$\int_0^\infty K_\nu(ax)I_\nu(bx)\cos(cx)dx=\frac{1}{2\sqrt{ab}}Q_{\nu-1/2}(\frac{a^2+b^2+c^2}{2ab})$,$\qquad{\rm Re}(a)>|{\rm Re}(b)|$, $c>0$, ${\rm Re}(\nu)>-1/2$

I already know the result of this integral because it is in Gradshtein & Ryzhik's book. Sometimes the integration is with respect to the order of the Bessel function.

Both Mathematica 8 and Maple 15 cannot do this kind of integrals. When the integral involves two Bessel functions or two other special functions, Mathematica and Maple usually cannot do even if the integral has a closed-form result. My questions are as follows:

1. Is there any general theory about how to do this kind of integrals? (I wonder how the authors of that book did the above integral.)

2. I know there is a Mathematica package "HolonomicFunctions." It seems that this package can help, but it does not seem very straightforward to obtain the final result. This package can verify the integrals with already known results, but can it do new integrals as above? Are there any better ways to deal with these integrals by computer?