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](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?



How to do integrals involving two Bessel functions and another function?
The best tool for trying to deal with such integrals is Fredéric Chyzak's [MGfun package](http://algo.inria.fr/chyzak/mgfun.html) (available as part of the [Algolib](http://algo.inria.fr/libraries/) library).


For your example, you should get a system of differential equations (for the integrand) for $a,b,c$ and $x$; you can leave $\nu$ as a parameter, or get a (mixed) difference equation for it. Then using this package, you can try to do elimination, which will give you a new system for the answer. 


Note that this method is a vast generalization of using the Meijer G-function as an intermediary (since MeijerG is the most general function whose series expansion / asymptotic expansion at 0 has coefficients satisfy a recurrence of order **1**). 


