I have been trying to evaluate the following integral involving modified Bessel functions:
\begin{align} \int r I_1 (r) K_1(r) dr \end{align}
This integral has an explicit expression given in the Wolfram documentation in terms of other Bessel functions. Granted, this expression has divergent components for integer-order Bessel functions but still has a well-defined limit at least for the above case. After performing these manually, one gets that the integral evaluates to
1/2 (-1 + r BesselI[0, r] (r BesselK[0, r] + BesselK[1, r]) +
r BesselI[1, r] (-BesselK[0, r] + r BesselK[1, r]))
However, when evaluating in Mathematica, this returns
In[1]:= Integrate[r BesselI[1, r] BesselK[1, r], r]
Out[1]= MeijerG[{{1, 3/2}, {}}, {{1, 2}, {0, 0}}, r, 1/2]/(4 Sqrt[\[Pi]])
Furthermore, this Meijer-G function refuses to simplify with FunctionExpand.
Is there a way get Mathematica to simplify the output of the integration? As far as I can see, the Mathematica internals should already know how to perform this simplification (or at least return the indefinite integral not in Meijer-G form).
Edit:
Furthermore, if one tries to evaluate
Integrate[z BesselI[\[Nu], a z] BesselK[\[Nu], a z], z]
as in the docs, Mathematica doesn't actually apply the rule discussed above, instead expressing in terms of hypergeometricPFQ functions. Why isn't this integration rule ever being used?
Series[MeijerG[{{1, 3/2}, {}}, {{1, 2}, {0, 0}}, r, 1/2]/( 4 Sqrt[\[Pi]]), {r, 0, 2}]
produces two "Power::infy: Infinite expression 1/0 encountered." and $$ \frac{G_{2,4}^{2,2}\left(0,\frac{1}{2}| \begin{array}{c} 1,\frac{3}{2} \\ 1,2,0,0 \\ \end{array} \right)}{4 \sqrt{\pi }}+\text{ComplexInfinity} r+\text{ComplexInfinity} r^2+O\left(r^3\right).$$ $\endgroup$