Find the first term expansion of the solutions of the following problem that is valid for large $t$, $\epsilon y'' +y' +y =0; \ \ \ t>0$
Regular expansion: Let us assume \begin{eqnarray} y(t) \sim y_0(t) + \epsilon y_1(t)+ ... . ...(1) \end{eqnarray} From the given equation we have \begin{eqnarray} \label{equ: substitution outer_expansion} y(t) \sim \epsilon (y_0''+ \epsilon y_1'' +...)+ (y_0'+ \epsilon y_1' +...)+(y_0+ \epsilon y_1 +...) \end{eqnarray} Then, \begin{eqnarray*} O(1) : \begin{cases} y_0'+y_0=0\\ y_0(0)=0 \end{cases} \end{eqnarray*}
Solving we get, $y_0(t)=c_1 e^{-t}$, where $c_1$ is a arbitrary constant.Using the initial condition we get, \begin{equation*} y_0(t)=0 \end{equation*}
Also, \begin{eqnarray*} O(\epsilon) : \begin{cases} y_0''+y'_1+y_1=0 \implies y'_1+y_1=0\\ y_1(0)=0 \end{cases} \end{eqnarray*}
Similarly we get \begin{equation*} y_1(t)=0 \end{equation*} Hence from (1) we get our solution, \begin{equation*} y(t) \sim 0 \end{equation*} Exact Solution: \begin{equation*} y(t) = \frac{\epsilon e^{\dfrac{(\sqrt{1-4 \epsilon}+1)x}{2\epsilon} } (e^{\dfrac{(\sqrt{1-4\epsilon})x}{\epsilon}}-1) }{\sqrt{1-4 \epsilon}} \end{equation*} Multiple-Scale Expansion: Consider two time scales $t_1=t$ and $t_2=\epsilon ^\alpha t$. \begin{equation*} \dfrac{d}{dt} \to \dfrac{dt_1}{dt} \dfrac{\partial }{\partial t_1} +\dfrac{dt_2}{dt} \dfrac{\partial }{\partial t_2}= \dfrac{\partial }{\partial t_1} + \epsilon^\alpha \dfrac{\partial }{\partial t_2} \end{equation*} Substituting this into the given equation we get, \begin{equation} \epsilon \bigg (\partial^2_{t_1} + 2 \epsilon^\alpha \partial_{t_1} \partial_{t_2}+ \epsilon^{2\alpha} \partial^2_{t_2}\bigg )y+ \bigg( \partial_{t_1}+ \epsilon^\alpha \partial_{t_2} \bigg )y+y=0 ...(2) \end{equation} where \begin{equation} y=0 \ \text{and} \ \bigg( \partial_{t_1}+ \epsilon^\alpha \partial_{t_2} \bigg )y=1, \ \text{for} \ t_1=t_2=0. \end{equation} Consider the power series expansion of the form
\begin{eqnarray} y \sim y_0(t_1, t_2) + \epsilon y_1(t_1, t_2)+ ... \end{eqnarray}
Substituting this into (2) yields the following
\begin{equation} \epsilon \bigg (\partial^2_{t_1} + 2 \epsilon^\alpha \partial_{t_1} \partial_{t_2}+ \epsilon^{2\alpha} \partial^2_{t_2}\bigg )(y_0+ \epsilon y_1 +...)+ \bigg( \partial_{t_1}+ \epsilon^\alpha \partial_{t_2} \bigg )(y_0+ \epsilon y_1 +...)+(y_0+ \epsilon y_1 +...)=0 \end{equation} \begin{eqnarray*} O(1) : \begin{cases} ( \partial_{t_1}+1)y_0=0,\\ y_0=0, \ \partial_{t_1}y_0=1 \ \ \text{at} \ \ t_1=t_2=0 \end{cases} \end{eqnarray*} The general solution of the problem is \begin{eqnarray} y_0= c(t_2) e^{t_1} \end{eqnarray}
Using the initial condition $y_0(0,0)=0$, we get $c(0)=0$. \ Balancing we get $\alpha=1$. \begin{eqnarray*} O(\epsilon) : \begin{cases} ( \partial_{t_1}+1)y_1=-(\partial^2_{t_1}+\partial_{t_2})y_0,\\ y_1=0, \ \partial_{t_1}y_1=- \partial_{t_2}y_0 \ \ \text{at} \ \ t_1=t_2=0 \end{cases} \end{eqnarray*} Then we get \begin{equation} ( \partial_{t_1}+1)y_1=-(c(t_2) +c'(t_2))e^{t_1} \end{equation} whose general solution is ,
\begin{equation} y_1= \bigg( (-c(t_2)-c'(t_2))t_1+p(t_2) \bigg)e^{-t_1} \end{equation} Using $y_1(0,0)=0$ ,we get $ p(0)=0 $. I am stuck here. Is there any secular term. If there is any secular term, how to prevent this.
First, I am assuming that $y(0)=0,y'(0)=1$. This isn't very clearly spelled out in your question, but what you wrote only really makes any sense under this assumption.
That said, your error is in how you selected $\alpha$ in your multiple time scales calculation. In actuality, the first time scale that needs to be resolved is the fast one, $\epsilon^{-1} t$ and then after that $t$. That is, there is a time scale faster than $t$ in the problem, which means that you really wanted to have $t_1=\epsilon^\alpha t,t_2=\epsilon^\beta t$. The equation becomes
$$\epsilon \left ( \epsilon^{2\alpha} \partial^2_{t_1} + 2 \epsilon^{\alpha+\beta} \partial_{t_1} \partial_{t_2} + \epsilon^{2\beta} \partial^2_{t_2} \right ) y + \left ( \epsilon^\alpha \partial_{t_1} + \epsilon^\beta \partial_{t_2} \right ) y + \epsilon^0 y =0.$$
Now there are numerous exponents on $\epsilon$ floating around: $1+2\alpha,1+\alpha+\beta,1+2\beta,\alpha,\beta$ and $0$. We can decide that $\alpha<\beta$, so that $t_1$ is the faster of the two scales. So you also have $1+2\alpha<1+\alpha+\beta<1+2\beta$. Therefore the exponents corresponding to the fastest scale have to be some two of $1+2\alpha,\alpha$ and $0$ (it cannot be all three at once, since $1+2(0) \neq 0$).
There are three such pairs: $\{ 1+2\alpha,\alpha \},\{ 1+2\alpha,0 \}$ and $\{ \alpha,0 \}$. For consistency, one must choose two of them to be equal such that those two are less than the third one.
Cases:
Plugging in $\alpha=-1$ you obtain
$$\partial^2_{t_1} y + \partial_{t_1} y = 0$$
which you can solve. Then the remaining exponents are $\beta,2\beta+1$ and $0$, and the correct balance becomes $\beta=0$. At this point, two uncompensated terms with exponent $1$ are left behind, which is a good sign. You get
$$2\partial_{t_1} \partial_{t_2} y + \partial_{t_2} y + y = 0.$$
You can use these two equations and the initial conditions to determine the leading order solution, you get a combination of a neither-fast-nor-slow decaying exponential and a fast decaying exponential.
An alternative way of getting a very similar result is boundary layer theory. Near the initial point, the $y''$ term must contribute, unless $y(0)+y'(0)=0$ (there are not enough constants of integration otherwise). So you zoom in there by considering $t=\epsilon s$ (the right scale so that $y''$ can balance with $y'$) and get an "inner" solution, controlled by $z''+z'=0$. You assume that away from the initial point, the $y''$ no longer contributes, so that the outer solution has $x'+x=0$. The composite of the two is given by matching the inner and outer solutions, adding them together, and then "subtracting the overlap". In this case the matching condition is $\lim_{s \to \infty} z(s)=\lim_{t \to 0} x(t)=:L$, which determines the constant of integration for $x$, and then $y(t)$ is approximated by $z(t/\epsilon)+x(t)-L$.
Notably the boundary layer approximation is already not quite as good as the one we got above, because it does not take into account that the fast scale is slightly slower than the inner solution suggests (which is a result of the $y$ term providing some "inertia"). I think this can be improved by applying regular perturbation theory to the inner and outer solutions but it sounds like a lot more effort than multiple scales analysis or the WKB method.