I want to find the first two leading terms of the expansion of $\int_0^{2\pi}ae^{x\cos a}da$
Well in $[0,2\pi]$ $\cos a$ has has maxima $0,2\pi$ so I rewrite the integral to $\int_0^{\epsilon}ae^{x\cos a} da+\int_{2\pi-\epsilon}^{2\pi}ae^{x\cos a}da$
$\cos a=1-\frac{a^2}{2}+O(a^4)$
How can I continue using Laplace Method?
You are on the right track. Essentially we need to approximate the integral near each maximum of $\cos a$. As you noted, there are two, one at $a = 0$ and $a = 2\pi$. We then fix $\epsilon > 0$ small and write the integral as
$$ \int_0^{2\pi} a e^{x\cos a}\,da = \int_0^\epsilon a e^{x\cos a}\,da + \int_{2\pi-\epsilon}^{2\pi} a e^{x\cos a}\,da + \int_{\epsilon}^{2\pi-\epsilon} a e^{x\cos a}\,da. \tag{1} $$
Our future calculations will verify that last integral here will not contribute to the asymptotic expansion as $x \to \infty$. For now we make note of the estimate
$$ \begin{align} \left|\int_{\epsilon}^{2\pi-\epsilon} a e^{x\cos a}\,da\right| &\leq \int_{\epsilon}^{2\pi-\epsilon} a e^{x\cos \epsilon}\,da \\ &= (2\pi^2-2\pi\epsilon) e^{x\cos \epsilon}. \tag{2} \end{align} $$
So, after replacing $a$ by $2\pi-a$ in the second integral in $(1)$, we know that
$$ \int_0^{2\pi} a e^{x\cos a}\,da = 2\pi \int_0^\epsilon e^{x\cos a}\,da + O\left(e^{x\cos\epsilon}\right) \tag{3} $$
as $x \to \infty$.
Note: We didn't have to make this change of variables to combine the two integrals here. We could have applied the next steps to each of the integrals in $(1)$ individually then summed the results.
Now, over the interval $[0,\epsilon]$ the integrand on the right-hand side of $(3)$ only has a maximum at $a = 0$, and near there we have
$$ \cos a = 1 - \frac{a^2}{2} + O(a^4). $$
This suggests the change of variables
$$ \cos a = 1 - b, $$
which yields
$$ \int_0^\epsilon e^{x\cos a}\,da = \frac{e^x}{\sqrt{2}} \int_0^{1-\cos\epsilon} b^{-1/2} \left(1-\frac{b}{2}\right)^{-1/2} e^{-xb}\,db. $$
The binomial theorem tells us that
$$ \left(1-\frac{b}{2}\right)^{-1/2} = \sum_{n=0}^{\infty} \binom{-1/2}{n}\left(-\frac{b}{2}\right)^n, $$
and we may therefore obtain the asymptotic expansion of the integral by appealing to Watson's lemma:
$$ \begin{align} &\frac{e^x}{\sqrt{2}} \int_0^{1-\cos\epsilon} b^{-1/2} \left(1-\frac{b}{2}\right)^{-1/2} e^{-xb}\,db \\ &\qquad \approx \frac{e^x}{\sqrt{2}} \sum_{n=0}^{\infty} \binom{-1/2}{n} \left(-\frac{1}{2}\right)^n \Gamma\left(n+\frac{1}{2}\right) x^{-n-1/2} \end{align} $$
as $x \to \infty$.
The error term in $(3)$, namely $O(e^{x\cos\epsilon})$, is dominated by each term in the above asymptotic series as $x \to \infty$ since they are all of the form $e^x x^\lambda$. It follows that the error does not contribute at all to the asymptotic series. We then conclude that
$$ \int_0^{2\pi} a e^{x\cos a}\,da \approx \sqrt{2}\pi e^x \sum_{n=0}^{\infty} \binom{-1/2}{n} \left(-\frac{1}{2}\right)^n \Gamma\left(n+\frac{1}{2}\right) x^{-n-1/2} $$
as $x \to \infty$.
