How to solve this set of definite integrals? I'm looking to solve the integrals
$$
I_n=\int_0^\infty \frac{x^ne^{-2ax}}{\sqrt{x^4+1}}dx.
$$
for $a>0$.
This should be reduced to finding $I_0$ as numerically I find that
$$
I_n=\left(-2\right)^{-n}\frac{d^nI_0}{da^n}.
$$
Mathematica gives a result in terms of the MeijerG function $G^{51}_{15}\left(x|^{\frac{1}{4}}_{-\frac{1}{2},-\frac{1}{4},-\frac{1}{4},0,\frac{1}{4}}\right)$ but I cannot manage to prove it. I have looked for integral representations of this and hypergeometric functions in Gradshteyn but couldn't find anything useful.
Any help is much welcomed.
Thanks

One way to get the result is to use the Mellin transform of the integral considered as a function of $a$:
\begin{align}
 I_0&=\int_0^\infty \frac{e^{-2ax}}{\sqrt{x^4+1}}\,dx\\
 \mathcal{M}\left[ I_0\right]&=\int_0^\infty a^{s-1}\,da\int_0^\infty \frac{e^{-2ax}}{\sqrt{x^4+1}}\,dx\\
 &=\int_0^\infty \frac{1}{\sqrt{x^4+1}}\,dx \int_0^\infty a^{s-1}e^{-2ax}\,da\\
 &=2^{-s}\Gamma(s)\int_0^\infty \frac{x^{-s}}{\sqrt{x^4+1}}\,dx
\end{align}
which is valid for $\Re(s)>0$. Now, changing $x=y^{1/4}$ in the integral, one obtains
\begin{align}
 \mathcal{M}\left[ I_n\right]&=2^{-s-2}\Gamma\left(s\right)\int_0^\infty\frac{y^{\frac{-s-3}{4}}}{\sqrt{y+1}}\,dy\\
 &=2^{-s-2}\Gamma(s)B\left(\frac{-s+1}{4},\frac{s+1}{4}\right)\\
 &=\frac{2^{-s}}{4\sqrt{\pi}}\Gamma\left(s\right)\Gamma\left(\frac{-s+1}{4}\right)\Gamma\left(\frac{s+1}{4}\right)
\end{align}
(The integral is the Mellin transform of $1/\sqrt{y+1}$ taken at $(n-s+1)/4$). This result is valid for $0<\Re(s)<1$. Taking the inverse transform, 
\begin{equation}
 I_0=\frac{1}{4\sqrt{\pi}}\frac{1}{2i\pi}\int_{\sigma-i\infty}^{\sigma+i\infty}(2a)^{-s}\Gamma\left(s\right)\Gamma\left(\frac{-s+1}{4}\right)\Gamma\left(\frac{s+1}{4}\right)\,ds
\end{equation}
with $0<\sigma<1$. Changing $s=4t$, one can express
\begin{equation}
 I_n=\frac{1}{\sqrt{\pi}}\frac{1}{2i\pi}\int_{\sigma-i\infty}^{\sigma+i\infty}(16a^4)^{-t}\Gamma\left(4t\right)\Gamma\left(t+\frac{1}{4}\right)\Gamma\left(\frac{1}{4}-t\right)\,dt
\end{equation} 
Expanding $\Gamma(4t)$ using the duplication formula:
\begin{equation}
 I_0=\frac{\sqrt{2}}{8\pi^2}\frac{1}{2i\pi}\int_{\sigma-i\infty}^{\sigma+i\infty}(\frac{a^4}{16})^{-t}\Gamma\left(t\right)\Gamma^2\left( t+\frac{1}{4} \right)\Gamma\left( t+\frac{1}{2} \right)\Gamma\left( t+\frac{3}{4} \right)\Gamma\left(\frac{1}{4}-t\right)\,dt
\end{equation} 
With $b_1=\frac{1}{4},a_1=1,a_2=\frac{3}{4},a_3=\frac{3}{4},a_4=\frac{1}{2},a_5=\frac{1}{4}$ and $\sigma=1/8$, we can express the result using the integral representation of the Meijer function DLMF:
\begin{equation}
 I_0=\frac{\sqrt{2}}{8\pi^2}{G^{1,5}_{5,1}}\left(\left.\frac{16}{a^4}\right|{1,\frac{3}{4},\frac{3}{4},\frac{1}{2},\frac{1}{4}\atop \frac{1}{4}}\right)
\end{equation} 
which, using these identities, can be transformed into
\begin{align}
 I_0&=\frac{\sqrt{2}}{8\pi^2}{G^{5,1}_{1,5}}\left(\left.\frac{a^4}{16}\right|{\frac{3}{4}\atop 0,\frac{1}{4},\frac{1}{4},\frac{1}{2},\frac{3}{4}}\right)\\
 &=\frac{\sqrt{2}}{32\pi^2}a^2{G^{5,1}_{1,5}}\left(\left.\frac{a^4}{16}\right|{\frac{1}{4}\atop -\frac{1}{2},-\frac{1}{4},-\frac{1}{4},0,\frac{1}{4}}\right)
\end{align} 
