Solve this differential equation using specific method! Consider the following problem:
$$\frac{\text{d}^2y}{\text{d}t^2} + \epsilon\frac{\text{d}y}{\text{d}t} + 1= 0\ \ \ \ \text{for}\ \  t \geq 0, y(0) = 0, y'(0)=0$$
a) Compute an approximate solution by substitution method:
 $y(t) = y_0(t) + \epsilon\ y_1(t) + \epsilon^2y_2(t) + \ldots \text{up to the order}\  \epsilon^2$;
b) Compute the exact solution
C) Using the Taylor's expansion of $e^-t$ at $t=0$, check that what you have seen in a) can be seen from the exact solution in b)
I know how to find an exact solution, but I have hard time using such substitution method to approximate the exact solution. Any help?

Substitute the series form of $y$ into the ODE, collect the $O(1)$, $O(\epsilon)$, and $O(\epsilon^2)$ terms, and solve each of these problems (ODEs) in this order. (Note, if  you are interested in perturbation theory of differential equations then Carl Bender's Mathematical Physics courses on Youtube is highly recommended.)
ADDED FULL SOLUTION:
$O(1)$ (unperturbed) problem:
$$
y_0''+1=0,
$$
which is solved to $y_0=-t^2/2+at+b$, but $a=b=0$ because we fit $y_0$ to the initial conditions.
$O(\epsilon)$ problem:
$$
\epsilon(y_1''+y_0')=0.
$$
Since $\epsilon\ne0$ and $y_0'=-t$ we can solve this easily for $y_1$ and get $y_1=t^3/6$. Note that the integration constants in $y_1$ are (again) zero because this time we fit the $O(\epsilon)$ solution, that is $y_0+\epsilon y_1$, to the initial conditions.
$O(\epsilon^2)$ problem:
$$
\epsilon^2(y_2''+y_1')=0.
$$
Along the same lines as above we get $y_1'=t^2/2$ and $y_2=-t^4/24$.
So putting everything together we finally have
$$
y = -\frac{t^2}{2}+\epsilon\frac{t^3}{6}-\epsilon^2\frac{t^4}{24} + O(\epsilon^3).
$$
The exact solution, which you have found too, as you claim, reads
$$
y = \frac{1}{\epsilon^2} - \frac{t}{\epsilon} - \frac{e^{-\epsilon t}}{\epsilon^2},
$$
and if you expand that in $\epsilon$ you will get the above series up to $O(\epsilon^2)$.
