A ﬁrst-term expansion of the solution of the ODE 
Get a ﬁrst-term expansion of the solution of $$εy''+y'+ y = 0, \text{where } y(0) = 0, y'(0) = 1$$ that is valid for large values of $t$.

I tried the usual approach with two time scales $t$ and $s=ε^a t$. Assuming $$y \sim y_0 + ε^b y_1$$ I get stuck immediately after writing it out since I can't see any obvious choice of $a$ and $b$.
Added: here is how I wrote it out: $$ε (\partial^2_t+2ε^a\partial_t \partial_s+ε^{2a}\partial^2_s)(y_0 + ε^b y_1) + (\partial_t + ε^a \partial_s)(y_0 + ε^b y_1) + (y_0 + ε^b y_1)=0$$
If we set $a=1$ and $b=1$, then $O(1)$ gives $$\partial_t y_0 + y_0 = 0$$ which can not satisfy our initial conditions.
So we want some other terms to be in $O(1)$ but both possible choices, namely $a = -1/2$ and $a = -1$, give us $1/ε$ in front of some terms.
 A: $$\epsilon y''+y'+y=0 \quad \begin{cases} y(0)=0\\y'(0)=1\end{cases}$$
NOTE :
In special case $\epsilon=0$ the solution of $y'+y=0$ is $y=ce^{-t}$ 
The condition $y(0)=0$ implies $c=0$ and the solution is trivial : $y(x)=0$
So, with the specified conditions, the solution of $\epsilon y''+y'+y=0$ tends to $y(x)=0$ when $\epsilon\to 0$. But the condition $y'(0)=1$ is ignored.
BRUT FORCE SOLVING METHOD :
With the eigen values $\epsilon r^2+r+1=0 \quad\to\quad r=\frac{-1\pm\sqrt{1-4\epsilon}}{2\epsilon}$ the general solution is :
$$y=c_1 e^{r_1t}+c_2 e^{r_2t} \quad 
\begin{cases}
r_1=\frac{-1+\sqrt{1-4\epsilon}}{2\epsilon} \sim -1-\epsilon\\
r_2=\frac{-1-\sqrt{1-4\epsilon}}{2\epsilon} \sim -\frac{1}{\epsilon}
\end{cases}$$
With conditions $y(0)=0=c_1+c_2$ and $y'(0)=1=r_1c_1+r_2c_2$ leading to $c_1=-c_2=\frac{\epsilon}{\sqrt{1-4\epsilon}}$
$$y=\frac{\epsilon}{\sqrt{1-4\epsilon}} \left( e^{\frac{-1+\sqrt{1-4\epsilon}}{2\epsilon} t}- e^{\frac{-1-\sqrt{1-4\epsilon}}{2\epsilon} t}\right)$$
$$y \sim \epsilon\left( e^{-t}- e^{-\frac{t}{\epsilon}}\right)$$
Bringing it back into the ODE and into the two conditions shows that it agrees for the first approximate.
This approch makes clear the difficulty when it is assumed that $y(t)\sim y_0(t)+\epsilon^b y_1(t)$ :
First, we have $y_0(t)=0$ .  Second, we have $\begin{cases} b=1\\ y_1(t)=e^{-t}\end{cases}$ . This is OK up to this point.
But the term $-\epsilon e^{-\frac{t}{\epsilon}}$is missing. This is acceptable for $\epsilon>0$ and $t$ large because if $\frac{t}{\epsilon}$ is large, the missing term $\epsilon e^{-\frac{t}{\epsilon}} \sim 0$.
On the other hand, keeping only the term $\epsilon e^{-t}$ doesn't satisfy the condition $y'(0)=1$.
Nevertheless, on the common sens of perturbation, I think that an acceptable answer to the question is :
$$y(t)\sim\epsilon e^{-t}$$
