Can we find the asymptotic behavior of this $f(x) =\int_{0}^{\infty}\frac{u^2}{1+\frac{e^{u^2}}{x}}du$?

In a first step, we perform the substitution $v=e^{u^2}$. we obtain $$ f(x) =\int_{0}^{\infty}\frac{u^2}{1+\frac{e^{u^2}}{x}}du = \int_1^\infty \frac{(\ln v)^{1/2}}{2 v(1+v/x)} dv \;.$$

Next, we take a look at the function $$ f_0(x) = \int_1^x \frac{(\ln v)^{1/2}}{2 v } dv =\frac13 (\ln x)^{3/2}$$ The intuition for this is that the integral for $f(x)$ is dominated for $v \lesssim x$. In this region the integral is approximately given by the expression above.

Indeed, we have $$f(x) - f_0(x) = -\int_1^x \frac{(\ln v)^{1/2}}{2(x+v)} dv+ \int_x^\infty \frac{(\ln v)^{1/2}}{2 v(1+v/x)} dv;$$ with the estimates $$\int_1^x \frac{(\ln v)^{1/2}}{2(x+v)} dv < \int_1^x \frac{(\ln x)^{1/2}}{2x} dv <\frac{(\ln x)^{1/2}}{2} $$ and $$\int_x^\infty \frac{(\ln v)^{1/2}}{2 v(1+v/x)} dv < \frac{x}{2}\int_x^\infty \frac{(\ln v)^{1/2}}{ v^2} dv = \frac12 (\ln x)^{1/2} + \frac{x}{4} \int_1^\infty \frac{1}{v^2(\ln v)^{1/2}} dv \\< \frac12 (\ln x)^{1/2} +\frac{x}{4 (\ln x)^{1/2}} \int_1^\infty \frac{1}{v^2} dv < \frac12 (\ln x)^{1/2} + \frac{1}{4 (\ln x)^{1/2}}\;. $$

With this, we have established that $$f(x) = \frac{1}{3} (\ln x)^{3/2} + O(\ln^{1/2} x).$$

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). $$

We have a linear upper bound: $$f=\int_0^\infty\frac{xu^2e^{-u^2}}{1+xe^{-u^2}}du\le x\int_0^\infty u^2e^{-u^2} du.$$But we also have a $O(1)$ lower bound: $$f\ge\frac{x}{1+xe^{-1}}\int_0^1u^2e^{-u^2}du=\frac{\int_0^1u^2e^{-u^2}du}{1/x+e^{-1}}.$$