Solving $y''(x) + \epsilon y'(x) + 1 = 0$ using power series We are given
\begin{align}
\begin{cases}
y''(x) + \epsilon y'(x) + 1 =0, \ 0 < \epsilon <<1\\
y(0)=0, \ y'(0)=1
\end{cases}
\end{align}
and asked to solve this using the solution form $y(x) = \sum_{n=0}^{\infty}\epsilon^{n}y_n(x)$.
Doing the known method for power series ODEs, we have
\begin{align}
&y(x) = \sum_{n=0}^{\infty}y_n(x) \epsilon^n\\
&y'(x) = \sum_{n=1}^{\infty}ny_n(x) \epsilon^{n-1} = \sum_{n=0}^{\infty}(n+1)y_{n+1}(x)\epsilon^{n}\\
&y''(x) = \sum_{n=2}^{\infty}n(n-1)y_n(x) \epsilon^{n-2} = \sum_{n=0}^{\infty}(n+2)(n+1)y_{n+2}(x)\epsilon^{n}.
\end{align}
So we must plug them into the system and find the general type $y_n$. I know how to do this in general. However, the present systems grinds to a halt, at least to my eyes.
Doing the substitution
\begin{align}
&\sum_{n=0}^{\infty} (n+2)(n+1)y_{n+2}(x)\epsilon^n + \sum_{n=0}^{\infty}(n+1)y_{n+1}(x) \epsilon^{n+1} + 1 =0\\
&1+ \sum_{n=0}^{\infty}(n+1)(n+1)y_{n+2}(x) \epsilon^{n} + \sum_{n=1}^{\infty}ny_n(x)\epsilon^{n} =0\\
&1 + 2y_2(x) + \sum_{n=1}^{\infty}\left\{ (n+2)(n+1)y_{n+2}(x) + ny_{n}(x) \right\}\epsilon^{n} = 0
\end{align}
and this means that
\begin{align}
\begin{cases}
y_2 = -\dfrac{1}{2}\\
y_{n+2} = \dfrac{-n}{(n+2)(n+1)}y_n
\end{cases}
\end{align}
but I think that leads to an algebraic fault, since by plugging $n=0$ at the second one we get $y_2(x) = 0$, but we found out that $y_2(x) = -\dfrac{1}{2}$.
Any thoughts on how to proceed?
EDIT: The answer is hinted to be
\begin{align}
y_{n}(x) = (-1)^{n} \left[ \dfrac{x^{n+1}}{(n+1)!} - \dfrac{x^{n+2}}{(n+2)!} \right]
\end{align}

One can solve directly
$$
y''+(ϵy'+1)=0\implies ϵy'(x)+1=(ϵy'(0)+1)e^{-ϵx}=(1+ϵ)e^{-ϵx}
\\
\implies
ϵy(x)+x=(1+ϵ)\frac{1-e^{-ϵx}}{ϵ}=(1+ϵ)\left(x-ϵ\frac{x^2}{2}+ϵ^2\frac{x^3}{3!}-ϵ^3\frac{x^4}{4!}+\dots\right)
\\
y(x)=x-\frac{x^2}{2}-ϵ\left(\frac{x^2}{2}-\frac{x^3}{3!}\right)+ϵ^2\left(\frac{x^3}{3!}-\frac{x^4}{4!}\right)+\dots
$$
So indeed, the reference solution
$$y(x)=\sum_{n=0}^\infty ϵ^ny_n(x),~~~y_n(x)=(-1)^n\left(\frac{x^{n+1}}{(n+1)!}-\frac{x^{n+2}}{(n+2)!}\right)
$$
is correct (contrary to my comment).

The x-derivatives of the perturbation series are, under the assumption of uniform convergence of all related series,
$$
y'(x)=\sum_{n=0}^\infty ϵ^ny_n'(x),~~~ y''(x)=\sum_{n=0}^\infty ϵ^ny_n''(x)
$$
After comparing coefficients, the resulting system is
$$
y_0''(x)+1=0,~~y_0(0)=0,~y_0'(0)=1\\
y_n''(x)+y_{n-1}'(x)=0,~~y_n(0)=y_n'(0)=0,~~n>0\\
$$
which indeed gives the same functions as solutions.
