An Irresistible integral: $\int_0^\infty \frac{x^{2m}}{(ax^2+b)^n}\mathrm dx=\frac{\pi}{2a^mb^{n-m-1}\sqrt{ab}}\frac{(2m-1)!!(2n-2m-3)!!}{(2m-2)!!}$

You may substitute $ax^2 + b = b/u$ to obtain

\begin{align*} \int_{0}^{\infty} \frac{x^{2m}}{(ax^2+b)^n} \, \mathrm{d}x &= \int_{0}^{1} \frac{\left( \frac{b}{a}(\frac{1}{u}-1) \right)^m}{\left(\frac{b}{u}\right)^n} \, \frac{\sqrt{b/a}}{2u^{3/2}(1-u)^{1/2}} \mathrm{d}u \\ &= \frac{(b/a)^{m+1/2}}{2 b^n} \int_{0}^{1} u^{n-m-3/2}(1-u)^{m-1/2} \, \mathrm{d}u \\ &= \frac{(b/a)^{m+1/2}}{2 b^n} \mathrm{B}\left(n-m-\tfrac{1}{2}, m + \tfrac{1}{2}\right). \end{align*}

Hope the computation is correct as I hurried a bit when writing things down...

Here is an alternative method where I will make use of a so-called Schwinger parametrisation. Such a parametrisation is particularly suited to integrals where a polynomial raised to a power appears in the denominator of the integrand.

For a positive, continuous function $\beta (x)$ observe that for $p > 0$ $$\frac{1}{\beta^p (x)} = \frac{1}{\Gamma (p)} \int_0^\infty u^{p - 1} e^{-u \beta(x)} \, du.$$ It is this observation that is known as the Schwinger parametrisation.

I will assume $a,b > 0$ and $m,n \in \mathbb{N}$ such that $n > m$. Choosing $\beta (x) = ax^2 + b$ and $p = n$ we have $$\frac{1}{(ax^2 + b)^n} = \frac{1}{\Gamma (n)} \int_0^\infty u^{n - 1} e^{-u(ax^2 + b)} \,du.$$ The integral $I(m, n; a, b)$ can thus be rewritten as \begin{align} I(m,n; a,b) &= \frac{1}{\Gamma (n)} \int_0^\infty \int_0^\infty x^{2m} u^{n - 1} e^{-u(ax^2 + b)} \, du \, dx\\ &= \frac{1}{\Gamma (n)} \int_0^\infty u^{n - 1} e^{-ub} \int_0^\infty x^{2m} e^{-uax^2} \, dx \, du, \end{align} after the order of integration has been changed.

Enforcing a substitution of $x \mapsto \sqrt{\dfrac{x}{ua}}$ leads to \begin{align} I(m,n;a,b) &= \frac{1}{2 \Gamma (n) a^{m + 1/2}} \int_0^\infty u^{n - m - 3/2} e^{-ub} \int_0^\infty x^{m - 1/2} e^{-x} \, dx\\ &= \frac{\Gamma (m + 1/2)}{2 \Gamma (n) a^{m + 1/2}} \int_0^\infty u^{n - m -3/2} e^{-ub} \, du. \end{align} Next, enforcing a substitution of $u \mapsto u/b$ yields \begin{align} I(m,n;a,b) &= \frac{\Gamma (m + 1/2)}{2 \Gamma (n) b^{n - m - 1/2} a^{m + 1/2}} \int_0^\infty u^{n - m -3/2} e^{-u} \, du\\ &= \frac{\Gamma (m + 1/2) \Gamma (n - m - 1/2)}{2 \Gamma (n) b^{n - m - 1/2} a^{m + 1/2}}, \end{align} or in terms of the beta function $$I(m,n;a,b) = \frac{1}{2b^n} \left (\frac{b}{a} \right )^{m + \frac{1}{2}} \operatorname{B} \left (n - m - \frac{1}{2}, m + \frac{1}{2} \right ),$$ in agreement with the result given by @Sangchul Lee. Granted, the substitution Sangchul Lee uses is pretty slick, and gets one to this point a lot quicker than making use of a Schwinger parametrisation, but it at least shows you an alternative approach to reaching the same point.

Note that in terms of central binomial coefficients the result can be expressed as $$I(m,n; a, b) = \frac{n \pi}{(2n - 2m - 1) 2^{2n} a^m b^{n - m - 1} \sqrt{ab}} \frac{\binom{2m}{m} \binom{2n - 2m}{n - m}}{\binom{n}{m}}.$$