Well, there are two contributions because the exponent is zero at $t=0$ and $t=\pi/2$. Let's consider $t=0$ first. In the immediate neighborhood of $t=0$ (we'll get to what that means in a bit), the exponent behaves as $x t^3$ so that as $x\to\infty$, we have
$$I(x) \sim \int_0^{\infty} dt \, e^{-x t^3} = \frac{\Gamma\left ( \frac{4}{3}\right )}{x^{1/3}} = \frac{\Gamma\left ( \frac{1}{3}\right )}{3 x^{1/3}} \quad (x\to\infty)$$
What is the size of the neighborhood? Well, we want $0 \lt x t^3 \lt \epsilon$ for some small $\epsilon$, so we have $0 \lt t \lt (\epsilon/x)^{1/3}$.
We also have a contribution in a neighborhood near $t=\pi/2$; we Taylor expand and get that
$$t^3 \cos{t} = -\frac{\pi^3}{8} \left ( t-\frac{\pi}{2}\right ) + O\left [ \left ( t-\frac{\pi}{2}\right )^2\right ]$$
Really, we are only interested in positive values of this expansion, as the exponent is positive through the integration region. If we look at only an immediate neighborhood near $t=\pi/2$, but with $t \lt \pi/2$, then we may approximate the contribution to the integral there as
$$\int_0^{\infty} dy \, e^{-\pi^3 x y/8} = \frac{8}{\pi^3 x}$$
It doesn't make sense to simply add these two terms together and declare them the leading behavior of $I(x)$ until we investigate the next leading behavior of the contribution at $x=0$. Note that the next contribution in the exponent is $x t^5/2$; within the interval of interest, this is $O(x^{-2/3})$, so we may Taylor expand this exponential term separately. The result is the following integral for the next contribution at $t=0$:
$$\frac12 x\int_0^{\infty} dt \;t^5 \, e^{-x t^3} = \frac1{6 x}$$
Note that this is $O(1/x)$ as is the leading contribution from $t=\pi/2$, so we may add these. The stated result follows.
The substitution $t=x(1-y)^{1/3}$, and Watson's lemma, give $$I(x)=\frac{x}{3}e^{x^3}\int_0^1(1-y)^{-2/3}e^{-x^3 y}\,dy\asymp\frac{e^{x^3}}{3}\sum_{n=0}^{(\infty)}\frac{\Gamma(n+2/3)}{\Gamma(2/3)x^{3n+2}}=e^{x^3}\left(\frac{1}{3x^2}+\frac{2}{9x^5}+\ldots\right).$$
Best Answer
Using the saddle point method with $a=0$, $b=+\infty$, $t_0=\mathrm{e}^{-1}$, $z=\mathrm{e}^n$, $p(t)=t\log t$ and $q(t)=1$, we find \begin{multline*} \int_0^{ + \infty } {x^{ - x} \mathrm{e}^{nx} \mathrm{d}x} = \mathrm{e}^n \int_0^{ + \infty } {\mathrm{e}^{ - \mathrm{e}^n t\log t} \mathrm{d}t} \\ \sim \exp \left( {\tfrac{{n - 1}}{2} +\mathrm{e}^{n - 1} } \right)\sqrt {2\pi } \left( {1 + \sum\limits_{k = 1}^\infty {\sqrt {\frac{{2\mathrm{e}}}{\pi }} \Gamma \left( {k + \tfrac{1}{2}} \right)\frac{{b_{2k} }}{{\mathrm{e}^{nk} }}} } \right), \end{multline*} where the coefficients are given as complex residues: $$ b_{2k} = \frac{1}{2}\mathop{\operatorname{Res}}\limits_{t = \mathrm{e}^{ - 1} }\left[ {(t\log t + \mathrm{e}^{ - 1} )^{ - k - 1/2} } \right] = \frac{1}{2}\mathrm{e}^{k - 1/2} \mathop{\operatorname{Res}}\limits_{s = 0} \left[ {((1 + s)\log (1 + s) - s)^{ - k - 1/2} } \right]. $$ For example, $$ b_2 = - \frac{1}{{12}}\sqrt {\frac{\mathrm{e}}{2}} ,\quad b_4 = - \frac{{23\mathrm{e}}}{{864}}\sqrt {\frac{\mathrm{e}}{2}} , $$ whence $$ \int_0^{ + \infty } {x^{ - x} \mathrm{e}^{nx} \mathrm{d}x} \sim \exp \left( {\tfrac{{n - 1}}{2} + \mathrm{e}^{n - 1} } \right)\sqrt {2\pi } \left( {1 - \frac{1}{{24}}\frac{1}{{\mathrm{e}^{n - 1} }} - \frac{{23}}{{1152}}\frac{1}{{\mathrm{e}^{2(n - 1)} }} + \cdots } \right). $$