Behavior of a function defined by two infinite integrals around the singularity point Bounty ending tomorrow:
Consider the following function defined by two infinite integrals:
$$
F(\epsilon) = \int_0^\infty 
 \frac{q^3 (1+q)e^{-2q}}{q^4+\epsilon^4} \, \mathrm{d} q
 +  \int_0^\infty \frac{\, G_q \, e^{-q}}{q^2(q^4+\epsilon^4)} \, \mathrm{d} q 
$$
where
$$
G_q := G \left( \left[ \left[\frac{1}{2} \right] , [\,]\right] , 
  \left[ \left[ \frac{9}{2},\frac{5}{2} \right], \left[\frac{3}{2}\right] \right] ,\frac{q^2}{4}  \right) \, , 
$$
with $G$ being the Meijer G-function.
The goal is to study the behavior of $F$ around the singularity point, i.e. as $\epsilon\to 0$.
Numerically, it can clearly be observed that $F$ scales logarithmically with $\epsilon$. I am wondering whether this behavior can be shown analytically via a rigorous analysis, e.g. using perturbation techniques. Your help or hints are highly appreciated and are most welcome. 
Thanks
H
 A: A technique similar to this will work. The tail at infinity can be neglected, and we have
$$ I_1 = \int_0^\infty \frac {q^3 (1 + q)} {q^4 + \epsilon^4} e^{-2q} dq \sim
\int_0^1 \frac {q^3} {q^4 + \epsilon^4} dq \sim
-\ln \epsilon.$$
The expansion for the G-function at zero is found from the residues:
$$G_q = G_{1,3}^{2,1}\left( \frac {q^2} 4 \middle| 
 {\frac 1 2 \atop \frac 5 2, \frac 9 2, \frac 3 2} \right) \sim\\
-\operatorname{Res}_{y=5/2}
 \frac {\Gamma\left( \frac 1 2 + y \right) \Gamma\left( \frac 5 2 - y \right)
   \Gamma\left( \frac 9 2 - y \right)}
  {\Gamma \left( - \frac 1 2 + y \right)}
 \left( \frac {q^2} 4 \right)^y = \frac {q^5} {16}.$$
Then, in the same way as for $I_1$,
$$I_2 = \int_0^\infty \frac {G_q} {q^2 \left( q^4 + \epsilon^4 \right)}
 e^{-q} dq \sim
\int_0^1 \frac {q^3} {16(q^4 + \epsilon^4)} dq \sim
-\frac {\ln \epsilon} {16},\\
F(\epsilon) = I_1 + I_2 \sim -\frac {17 \ln \epsilon} {16}.$$
To obtain the next term, a regularization can be applied to $I_1$ to get
$$I_1 \sim -\ln \epsilon +
\int_0^1 \left( \frac {q^3 (1 + q)} {q^4 + \epsilon^4} e^{-2q} -
 \frac 1 q \right) \bigg\rvert_{\epsilon=0} dq +
\int_1^\infty \frac {q^3 (1 + q)} {q^4 + \epsilon^4} e^{-2q}
 \bigg\rvert_{\epsilon=0} dq,$$
and, repeating the same procedure for $I_2$, we find
$$F(\epsilon) \sim -\frac {17 \ln \epsilon} {16} -
 \gamma - \ln 2 + \frac 1 2 +\\
 \int_0^1 \left( q^{-6} e^{-q} G_q - \frac 1 {16 q} \right) dq +
 \int_1^\infty q^{-6} e^{-q} G_q dq.$$
