I would appreciate help in getting Mathematica to solve an indefinite integral that is closely related to two that it can solve, or at least to understand why it cannot. Consider an indefinite integral that includes the half-integer Macdonald function, also known as the “modified spherical Bessel function of the second kind” and the “reduced Bessel function,” multiplied by $x^{\frac{1}{2}-h} e^{-a x^2}$. If $n=2$ and $h=0$ or $2$, Mathematica 7 and 13 can solve the integral: $$\int e^{-a x^2} x^{\frac{1}{2}-0} K_{\frac{1}{2}+2}(x b) b^{\frac{1}{2}+2} dx = \sqrt{\frac{\pi }{2}} \left(-\frac{3 e^{-x (b+a x)}}{x}-\frac{\left(6 a-b^2\right) e^{\frac{b^2}{4 a}} \sqrt{\pi } \text{erf}\left(\frac{b+2 a x}{2 \sqrt{a}}\right)}{2 \sqrt{a}}\right) \quad (1)$$ and $$\int e^{-a x^2} x^{\frac{1}{2}-2} K_{\frac{1}{2}+2}(x b) b^{\frac{1}{2}+2} dx =$$ $$\sqrt{\frac{\pi }{2}} \left(e^{-b x-a x^2} \left(-\frac{1}{x^3}-\frac{b}{x^2}+\frac{2 a}{x}\right)+2 a^{3/2} e^{\frac{b^2}{4 a}} \sqrt{\pi } \text{erf}\left(\frac{b+2 a x}{2 \sqrt{a}}\right)\right) \quad, \quad (2)$$ respectively. But (now shifting to InputForm)

Integrate[   E^(- a (x^2)) x^(1/2 - h) BesselK[1/2 + n, x b] b^(1/2 + n), x] /. 
  n -> 2 /. h -> 1

yields the response $$b^{5/2} \sqrt{\frac{\pi }{2}} \int \frac{e^{-b x-a x^2} \left(1+\frac{3}{b^2 x^2}+\frac{3}{b x}\right)}{\sqrt{x} \sqrt{b x}} \, dx $$ that indicates a failure to solve the integral.

One gets the same response if one directly integrates the above series expression for the half-integer Macdonald function

b^2*Sqrt[Pi/2]*Integrate[(E^((-b)*x - a*x^2)*(1 + 3/(b^2*x^2) + 3/(b*x)))/x, x]

(or attempts to integrate the series term-by-term) or uses the Meijer-G function form of the half-integer Macdonald function:

Integrate[(x^2 b^2)^(1/4)/(Sqrt[x] Sqrt[b]) E^(-x b)/E^-Sqrt[x^2 b^2]  E^(- a (x^2)) x^(1/2 - h) 
    *1/2 MeijerG[{{}, {}}, {{1/2 (1/2 + n), 1/2 (-(1/2) - n)}, {}}, 1/4 x^2 b^2] * b^(1/2 + n), x] /. n -> 2 /. 
          h -> 1 /. {1/Sqrt[x^2 b^2] -> 1/(x b)}

where one needs the two factors of one on the first line, (x^2 b^2)^(1/4)/(Sqrt[x] Sqrt[b]) and E^(-x b)/E^-Sqrt[x^2 b^2], and the substitution /. {1/Sqrt[x^2 b^2] -> 1/(x b)} on the last line in order for Mathematica to yield the right-hand sides of equations (1) and (2) when $h$ is set to 0 and 2, resp.

So my question is, since Mathematica can numerically integrate these expressions when $h=0$, 1, or 2, and since it can solve for the indefinite integral when $h$ takes on the bracketing values $h=0$ and $h=2$, why will it not solve for the indefinite integral at the intermediate value $h=1$? Is there some way to get it to do so? (It also fails to do so for even values $h=4$ and larger but can for positive powers of x in the series given by $h=-2$.)

Perhaps I should instead be thankful that Mathematica can solve for the indefinite integral for the bracketing values $h=0$ and $h=2$ since equations (1) and (2) are not tabled in Gradshteyn and Ryzhik (Section 5.5) nor Prudnikov, Brychkov, and Marichev (Section 1.12.2) and ask how it manages to do so given that it cannot integrate the series term-by-term.

Thanks much, Jack