How to do integrals involving two Bessel functions and another function? 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:


*

*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.)

*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?
 A: The best tool for trying to deal with such integrals is Fredéric Chyzak's MGfun package (available as part of the Algolib 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).  
A: Take a look at "A Treatise on the Theory of Bessel Functions" by Watson. There is a long chapter on integrating Bessel functions over the infinite range $0-\infty$. 
In addition, I think a Mellin transform approach could very well get you what you want. The idea is to find Mellin transforms of the functions in the integrand (most often separated into two parts) and use Parseval's theorem to write the integral as a contour integral. Then one can often move the contour over the poles of the integrand and generate a series representation of the integral, which can sometimes be identified as some known special function. 
You could look at the paper http://www.risc.jku.at/publications/download/risc_3924/AlgorithmicMellinTransform.pdf which briefly describes the method and shows a computer algebra technique for getting the final result. There is also a great, simple book by Fikoris called "Mellin Transform Method for Integral Evaluation." 
Good luck,
Tom
