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.
In Bender & Orszag, at least in my edition (the first, 1978), this is problem 6.73 on p. 313. On p. 308, problem 6.19(a), you are asked how many terms one can get using integration by parts. (The answer is one - the second term blows up at the $t=0$ limit.) Your answer is correct for the leading term.
The idea here is to use steepest descents. I will outline here how to go about setting up the steepest descent contours and setting up the integrals used to derive the asymptotic expansions. Really, this case follows closely Example 2 of Sec. 6.6 in B & O.
Recall that the whole idea behind steepest descent is to find an approximation to the complex integral
$$\int_C dz \, h(z) e^{x \rho(z)} $$
as $x \to \infty$. We do this by deforming $C$ to a steepest descent contour $C'$ along which we may end up getting integrals to which we may apply Laplace's method, e.g., decaying exponentials. So, again, consider the integral of interest:
$$\int_0^{\pi/4} dt \, \cos{x t^2} \, \tan^2{t} = \operatorname{Re} \int_0^{\pi/4} dt \, e^{i x t^2} \, \tan^2{t}$$
Our contour $C$ is just a line segment over the real axis. What would be the steepest descent contour $C'$? It is that over which $\operatorname{Im}{\rho(z)}$ is constant. Letting $z=u+i v$, this would mean that we want a contour over which $i z^2 = -2 u v+i (u^2-v^2) $ is constant.
At $z=0$, $\operatorname{Im}{(i z^2)}=0$, so we can use a contour defined by $u=v$. ($u=-v$ defines a steepest ascent contour, which does not deliver an integral useful for Laplace's method.) Thus, the contour $C_1$ coming from $z=0$ is $z=e^{i \pi/4} t$.
We will discuss a contour $C_3$ coming from the other endpoint at $z=\pi/4$ in a minute. We deform the original integral into a closed contour $\gamma=C_1+C_2-C_3-C$ as follows:
The contour $C_2$ is the top piece and vanishes as we go farther along the contours $C_1$ and $C_3$. Thus, by Cauchy's theorem we have
$$\int_0^{\pi/4} dt \, e^{i x t^2} \tan^2{t} = \int_{C_1} dz \, e^{i x z^2} \tan^2{z} - \int_{C_3} dz \, e^{i x z^2} \tan^2{z} $$
Now,
$$\int_{C_1} dz \, e^{i x z^2} \tan^2{z} = e^{i \pi/4} \int_0^{\infty} dt \, e^{-x t^2} \tan^2{\left ( e^{i \pi/4} t \right )} $$
Now, at $z=\pi/4$, $\operatorname{Im}{(i z^2)}=\pi^2/16$. Thus, along $C_3$, $u^2-v^2=\pi^2/16$, or $u=\sqrt{v^2+\pi^2/16}$ and
$$\rho(z) = i z^2 = i \frac{\pi^2}{16} - 2 v \sqrt{v^2+\frac{\pi^2}{16}} $$
As in the book, let $s= 2 v \sqrt{v^2+\frac{\pi^2}{16}}$ and so on, and we find that
$$\int_{C_3} dz \, e^{i x z^2} \tan^2{z} = i e^{i \pi^2 x/16} \frac{2}{\pi} \int_0^{\infty} ds \, \left (1+i \frac{16 s}{\pi^2} \right )^{-1/2} \tan^2{\left (\frac{\pi^2}{16}+i s \right )^{1/2}} \, e^{-x s}$$
I hope you can see where this is going. We expand separately along each contour for large $x$. At $C_1$, i.e. $t=0$, it should be clear that
$$\begin{align} e^{i \pi/4} \int_0^{\infty} dt \, e^{-x t^2} \tan^2{\left ( e^{i \pi/4} t \right )} &= i e^{i \pi/4} \int_0^{\infty} dt \, e^{-x t^2} [t^2+O(t^3)]\\ &= -i e^{i \pi/4} \frac14 \sqrt{\pi} x^{-3/2} + O(x^{-2}) \end{align}$$
At $C_3$, i.e., $t=\pi/4$, we have
$$\begin{align} i e^{i \pi^2 x/16} \frac{2}{\pi} \int_0^{\infty} ds \, \left (1+i \frac{16 s}{\pi^2} \right )^{-1/2} \tan^2{\left (\frac{\pi^2}{16}+i s \right )^{1/2}} \, e^{-x s} &= i e^{i \pi^2 x/16} \frac{2}{\pi} \int_0^{\infty} ds \, \left [1+O(s)\right ] e^{-x s} \\ &= i e^{i \pi^2 x/16} \frac{2}{\pi x} + O(x^{-2}) \end{align}$$
Finally, subtracting the two contributions and taking the real part, we get the first two terms of the expansion:
$$\int_0^{\pi/4} dt \, \cos{x t^2} \, \tan^2{t} = \frac{2}{\pi x} \sin{\left (\frac{\pi^2}{16} x \right )} - \frac14 \sqrt{\frac{\pi}{2}} x^{-3/2} + O(x^{-2})$$
The nice part is that we can keep adding on terms if we want.
Here is a graph verifying the above asymptotic expansion:
The blue plot is the log-log error between a numerical value of the integral and the first term in the expansion. The red plot is the log-log-error with both terms. Note that the first plot has a supremum that has slope -1.5, and the red has one that has slope -2.0. This indicates that the expansion above is correct.
Best Answer
First approach. We have \begin{align*} b_n & = \frac{{{\rm e}^{ - n} }}{{\Gamma (n)}}\int_0^{ + \infty } {\frac{{\Gamma (x + n)}}{{\Gamma (x + 1)}}{\rm e}^{ - x} {\rm d}x} \\ & = \frac{{{\rm e}^{ - n} }}{{\Gamma (n)}}\int_0^{ + \infty } {\frac{1}{{\Gamma (x + 1)}}{\rm e}^{ - x} \left( {\int_0^{ + \infty } {s^{x + n - 1} {\rm e}^{ - s} {\rm d}s} } \right)\!{\rm d}x} \\ & = \frac{1}{{\Gamma (n)}}\int_0^{ + \infty } {\frac{1}{{\Gamma (x + 1)}}\left( {\int_0^{ + \infty } {t^{x + n - 1} {\rm e}^{ - {\rm e}t} {\rm d}t} } \right)\!{\rm d}x} \\ & = \frac{1}{{\Gamma (n)}}\int_0^{ + \infty } {t^{n - 1} {\rm e}^{ - {\rm e}t} \left( {\int_0^{ + \infty } {\frac{{t^x }}{{\Gamma (x + 1)}}{\rm d}x} } \right)\!{\rm d}t} . \end{align*} Employing Ramanujan's formula $$ \int_0^{ + \infty } {\frac{{t^x }}{{\Gamma (1 + x)}}{\rm d}x} = {\rm e}^t - \int_{ - \infty }^{ + \infty } {\frac{{{\rm e}^{ - t{\rm e}^y } }}{{y^2 + \pi ^2 }}{\rm d}y} , $$ yields the exact expression \begin{align*} b_n & = \frac{1}{{\Gamma (n)}}\int_0^{ + \infty } {t^{n - 1} {\rm e}^{ - ({\rm e} - 1)t} {\rm d}t} - \frac{1}{{\Gamma (n)}}\int_0^{ + \infty } {t^{n - 1} {\rm e}^{ - {\rm e}t} \int_{ - \infty }^{ + \infty } {\frac{{{\rm e}^{ - t{\rm e}^y } }}{{y^2 + \pi ^2 }}{\rm d}y}\, {\rm d}t} \\ & = \frac{1}{{({\rm e} - 1)^n }} - \int_{ - \infty }^{ + \infty } {\frac{1}{{({\rm e} + {\rm e}^y )^n }}\frac{1}{{y^2 + \pi ^2 }}{\rm d}y} . \end{align*} Since $$ \int_{ - \infty }^{ + \infty } {\frac{1}{{({\rm e} + {\rm e}^y )^n }}\frac{1}{{y^2 + \pi ^2 }}{\rm d}y} \le \frac{1}{{{\rm e}^n }}\int_{ - \infty }^{ + \infty } {\frac{{{\rm d}y}}{{y^2 + \pi ^2 }}} = \frac{1}{{{\rm e}^n }}, $$ we indeed have $$ b_n \sim \frac{1}{{({\rm e} - 1)^n }} $$ as $n\to +\infty$.
Second approach. Changing the order of summation and integration yields $$ \sum\limits_{n = 1}^\infty {b_n z^n } = z\int_0^{ + \infty } {\frac{{{\rm d}x}}{{({\rm e} - z)^{x + 1} }}} = \frac{z}{{({\rm e} - z)\log ({\rm e} - z)}} $$ for sufficiently small $z$. Now note that $$ \frac{z}{{({\rm e} - z)\log ({\rm e} - z)}} = \frac{{\rm e} - 1}{{({\rm e} - 1) - z}} + H(z) $$ where $H(z)$ is holomorphic in the disc $|z|<\mathrm{e}$. The first term may be expanded as $$ \frac{{\rm e} - 1}{{({\rm e} - 1) - z}} = \sum\limits_{n = 0}^\infty {\frac{1}{{({\rm e} - 1)^{n} }}z^n } . $$ On the other hand, the $n$th Maclaurin series coefficient of $H(z)$ is $\mathcal{O}((\mathrm{e}-\varepsilon)^{-n})$ by the Cauchy–Hadamard theorem for any $\varepsilon>0$ as $n\to+\infty$. Thus $$ b_n \sim \frac{1}{{({\rm e} - 1)^n }} $$ as $n\to +\infty$.