Asymptotic evaluation of $\int_0^{\pi/4}\cos(x t^2)\tan^2(t)dt$ In Bender-Orszag's  Advanced Mathematical Methods for Scientists and Engineers on page 313 we encounter the following integral 

$$
I(x)=\int_0^{\pi/4}\cos(x t^2)\tan^2(t)dt
$$

and it is asked to derive the first two terms of the asymptotic expansion of $I(x)$ as $x\rightarrow +\infty$ . 
Setting $x=i y$ and analytically continuing at the end of the calcultion I was able to show that

$$
I(x)\sim 2\frac{\sin(\frac{x\pi^2}{16})}{\pi x}+\mathcal{O}(x^{-2})
$$

using standard Laplace method. Despite the fact that this seems to fit comparsion with numerical data, I am not satisfied by this approach for three reasons:

1) It is quiet cumbersome (at least for me) to derive higher order terms with my method. Is there a fast way to do so?
2) The question can be found in relation with "Method of steepest descent" so I suppose one should attack it by this method. How can it be applied? I somehow fail here
3) I am not totally satisfied by my analytic continuation argument. How can this made rigourous?

I am grateful to anyone who can shed light on any of my questions!
 A: 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.
A: I know that this may not be the direction that OP want, but we can revive @Olivier Oloa's idea to give an elementary solution. Taking integration by parts, we get
\begin{align*}
I_0(x)
:\!\!&= \int_{0}^{\frac{\pi}{4}} \cos(xt^2)\tan^2 t \, dt \\
&= \frac{2}{\pi}\frac{\sin\big(\frac{\pi^2}{16}x\big)}{x} + \frac{1}{2x} \underbrace{\int_{0}^{\frac{\pi}{4}} \sin(xt^2) \bigg( \frac{\tan^2 t}{t^2} - \frac{2\tan t\sec^2 t}{t} \bigg) \, dt}_{=:I_1(x)}.
\end{align*}
Now we want to apply a similar technique to $I_1(x)$ to extract higher order terms, but this technique need a modification in the following way: Write
\begin{align*}
I_1(x)
&= \int_{0}^{\frac{\pi}{4}} \sin(xt^2) \bigg( 1 + \frac{\tan^2 t}{t^2} - \frac{2\tan t\sec^2 t}{t} \bigg) \, dt - \int_{0}^{\frac{\pi}{4}} \sin(xt^2) \, dt.
\end{align*}
This modification takes care of the singularity that would otherwise have popped up under integration by parts. So the same technique applies to the first integral above and shows that it is of order $\mathcal{O}(x^{-1})$. So we get
$$ I_1(x) = - \int_{0}^{\frac{\pi}{4}} \sin(xt^2) \, dt + \mathcal{O}(x^{-1}). $$
On the other hand, using the substitution $u = xt^2$ and assuming $x > 0$, we have
\begin{align*}
\int_{0}^{\frac{\pi}{4}} \sin(xt^2) \, dt
&= \frac{1}{2\sqrt{x}} \int_{0}^{\frac{\pi^2}{16}x} \frac{\sin u}{\sqrt{u}} \, du \\
&= \frac{1}{2\sqrt{x}} \int_{0}^{\infty} \frac{\sin u}{\sqrt{u}} \, du - \frac{1}{2\sqrt{x}} \int_{\frac{\pi^2}{16}x}^{\infty} \frac{\sin u}{\sqrt{u}} \, du \\
&= \frac{1}{2\sqrt{x}} \sqrt{\frac{\pi}{2}} - \frac{1}{2\sqrt{x}} \bigg( \underbrace{\frac{4}{\pi}\frac{\cos \big(\frac{\pi^2}{16}x\big)}{\sqrt{x}} + \frac{1}{2} \int_{\frac{\pi^2}{16}x}^{\infty} \frac{\cos u}{u^{3/2}} \, du}_{=\mathcal{O}(x^{-1/2})} \bigg).
\end{align*}
Combining altogether, we get
$$ I_0(x) = \frac{2}{\pi}\frac{\sin\big(\frac{\pi^2}{16}x\big)}{x} - \frac{1}{4}\sqrt{\frac{\pi}{2}}x^{-3/2}+ \mathcal{O}(x^{-2}). $$
A: User Ron Gordon has already given a good answer. Here we will try to extract, streamline & generalize his ideas.
I) The underlying complex integral is of the form
$$\begin{align}I(a)~:=~&\int_0^a  \! \mathrm{d}z~ \exp\left(ix z^n\right) f(z),\cr 
a~\geq~&0, \qquad n~\in~\mathbb{N},\qquad x~>~0  ,\end{align}\tag{1}$$
where $f$ is a meromorphic function. Rather than using the steepest descent/stationary phase/saddle-point approximation, below method contains the following ingredients:

*

*Deform the integration contour (in principle not necessarily in the direction of steepest descent) into a sector of the complex $z$-plane, where the integrand is exponentially suppressed. Keep in mind residues from crossing poles of $f$.


*Use the integral representation for the Euler gamma function rather than Gaussian integration.


*Up until now the method is in principle exact. Finally Taylor expansion to the needed order of precision.
The complex integral (1) becomes
$$I(a)~=~ J(0)-\exp\left(ixa^n\right)J(a), \tag{2}$$
where
$$\begin{align} J(a)~:=~&\int_0^{\infty \exp\left(\frac{i\pi}{2n}\right)} \! \mathrm{d}z~  \exp\left[ix\left\{(z+a)^n-a^n \right\}\right] f(z+a) \cr
~\stackrel{(4)}{=}~&\frac{i}{n} \int_0^{\infty} \! \mathrm{d}u~  \exp\left(-xu\right) \frac{f\left( \sqrt[n]{iu+a^n} \right)}{(iu+a^n)^{\frac{n-1}{n}}} .\end{align}\tag{3} $$
Here we have used the substitution
$$z~=~\sqrt[n]{iu+a^n}-a 
\qquad \Rightarrow \qquad 
iu~=~(z+a)^n-a^n ,\qquad u~\geq~0. \tag{4} $$
II) For the rest of this answer let us specialize to OP's example with $n=2$ and $f(z)=\tan^2z$.
The integral associated with the upper endpoint $a=\frac{\pi}{4}$ yields
$$\begin{align} J\left(a\!=\!\frac{\pi}{4}\right)
~\stackrel{(3)}{=}~&\frac{i}{2} \int_0^{\infty} \! \mathrm{d}u~  \exp\left(-xu\right)  \left( \frac{4}{\pi} + \frac{32(\pi-1)}{\pi^3} iu + O(u^2) \right) \cr
~=~&\frac{2i}{\pi x} -\frac{16(\pi-1)}{x^2\pi^3} + O(x^{-3})  .\end{align}\tag{5} $$
The integral associated with the lower endpoint $a=0$ yields
$$\begin{align} J(a\!=\!0)~\stackrel{(3)}{=}~&
-\frac{1}{2} \exp\left(\frac{i\pi}{4}\right) \int_0^{\infty} \! \mathrm{d}u~  \exp\left(-xu\right) \frac{\tanh^2 \sqrt{u}}{\sqrt{u}}\cr
~=~&-\frac{1}{2} \exp\left(\frac{i\pi}{4}\right) \int_0^{\infty} \! \mathrm{d}u~  \exp\left(-xu\right) \left(u^{\frac{1}{2}} -\frac{2}{3} u^{\frac{3}{2}} + O(u^{\frac{5}{2}})\right)\cr
~=~&  -\frac{\sqrt{\pi}}{4 x^{\frac{3}{2}}}\exp\left(\frac{i\pi}{4}\right)\left( 1-x^{-1}  + O(x^{-2})\right) .\end{align}\tag{6} $$
This leads to the OP's sought-for expansion

$$ {\rm Re}I\left(a\!=\!\frac{\pi}{4}\right)~\stackrel{(2)+(5)+(6)}{=}~\frac{2}{\pi x} \sin{\left(\frac{\pi^2}{16} x \right)} 
- \frac{\sqrt{2\pi}}{8x^{\frac{3}{2}}} + O(x^{-2}) .\tag{7}$$

