Can we find the asymptotic behavior of this $f(x) =\int_{0}^{\infty}\frac{u^2}{1+\frac{e^{u^2}}{x}}du$? I encountered this function in Statistical Mechanics.  $$f(x) =\int_{0}^{\infty}\frac{u^2}{1+\frac{e^{u^2}}{x}}du$$ 
For $x=0$, we define its value to be zero.
I wanted to see it's asymptotic behavior in the limit x tending to $\infty$. Can we express the asymptotic behavior of this function, in the limit x tending to $\infty$, in terms of other known mathematical functions (if possible; in elementary functions)? For some physical reasons (which are irrelevant here), $x$ belongs to $[0,\infty)$.

1. (Not so illuminating) analytic expression. Assume for a moment that $0 < x < 1$. Then
\begin{align*}
f(x)
&= \int_{0}^{\infty} \frac{xu^2e^{-u^2}}{1 + xe^{-u^2}} \, du
 = \sum_{n=1}^{\infty} (-1)^{n-1} x^n \int_{0}^{\infty} u^2 e^{-nu^2} \, du \\
&= -\frac{\sqrt{\pi}}{4} \sum_{n=1}^{\infty} \frac{(-x)^n}{n^{3/2}}
 = -\frac{\sqrt{\pi}}{4} \operatorname{Li}_{3/2}(-x).
\end{align*}
The last function is analytic outside $(-\infty, -1]$, and hence this identity extends to all of $x \geq 0$ by the principle of analytic continuation. But this is not so useful when investigating the asymptotic bahavior of $f(x)$.
2. Asymptotic expansion. Write $\alpha = \log x$ and make the substitution $u = \sqrt{\alpha(v+1)}$. Then
\begin{align*}
f(x)
&= \frac{\alpha^{3/2}}{2} \int_{-1}^{\infty} \frac{\sqrt{1+v}}{1 + e^{\alpha v}} \, dv \\
&= \frac{\alpha^{3/2}}{2} \left( \int_{0}^{1} \sqrt{1-v} \, dv - \int_{0}^{1} \frac{\sqrt{1-v}}{1 + e^{\alpha v}} \, dv + \int_{0}^{\infty} \frac{\sqrt{1+v}}{1 + e^{\alpha v}} \, dv \right).
\end{align*}
This easily yields the following asymptotics
$$ f(x) = \frac{1}{3}(\log x)^{3/2} + \mathcal{O}\left( (\log x)^{1/2} \right). $$
For a better resolution, recall that the polylogarithm is defined as $\operatorname{Li}_s(z) = \sum_{n=1}^{\infty} \frac{z^n}{n^s}$ for $|z| < 1$. Then
$$ \frac{1}{1 + e^{\alpha v}} = -\operatorname{Li}_0(-e^{-\alpha v}), \qquad \frac{d}{dv} \operatorname{Li}_{s+1}(-e^{-\alpha v}) = - \alpha \operatorname{Li}_s(-e^{-\alpha v}) $$
and hence
\begin{align*}
\int_{0}^{\infty} \frac{\sqrt{1+v}}{1 + e^{\alpha v}} \, dv
&= -\int_{0}^{\infty} (1+v)^{1/2} \operatorname{Li}_0(-e^{-\alpha v}) \, dv \\
&= -\frac{\operatorname{Li}_1(-1)}{\alpha} - \frac{1}{2\alpha} \int_{0}^{\infty} \frac{\operatorname{Li}_1(-e^{-\alpha v})}{(1+v)^{1/2}} \, dv \\
&= -\frac{\operatorname{Li}_1(-1)}{\alpha} - \frac{\operatorname{Li}_2(-1)}{2\alpha^2} + \frac{1}{4\alpha^2} \int_{0}^{\infty} \frac{\operatorname{Li}_2(-e^{-\alpha v})}{(1+v)^{3/2}} \, dv
\end{align*}
and, in principle, the same argument can be applied to extract an asymptotic expansion up to any fixed order. Similarly,
\begin{align*}
\int_{0}^{1} \frac{\sqrt{1-v}}{1 + e^{\alpha v}} \, dv
&= -\int_{0}^{1} (1-v)^{1/2} \operatorname{Li}_0(-e^{-\alpha v}) \, dv \\
&= -\frac{\operatorname{Li}_1(-1)}{\alpha} + \frac{1}{2\alpha} \int_{0}^{1} \frac{\operatorname{Li}_1(-e^{-\alpha v})}{(1-v)^{1/2}} \, dv \\
&= -\frac{\operatorname{Li}_1(-1)}{\alpha} + \frac{\operatorname{Li}_2(-1) - \operatorname{Li}_2(e^{-\alpha})}{2\alpha^2} \\
&\qquad + \frac{1}{4\alpha^2} \int_{0}^{1} \frac{\operatorname{Li}_2(-e^{-\alpha v}) - \operatorname{Li}_2(-e^{-\alpha})}{(1-v)^{3/2}} \, dv
\end{align*}
and so on. Using the results above, we obtain a better asymptotics
$$ f(x) = \frac{1}{3} (\log x)^{3/2} + \frac{\pi^2}{24} \frac{1}{(\log x)^{1/2}} + \mathcal{O}\left( \frac{1}{(\log x)^{3/2}} \right). $$
