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}
 A: You methodology is basically correct but there are some slips in the algebra. I feel more comfortable writing $y_n(x)=a_nx^n.$ Then $y=a_0+\epsilon a_1x+\epsilon^2a_2x^2+...$ Then $y(0)=a_0=0$ and$y^{\prime}(0)=\epsilon a_1$, so $a_1=1/\epsilon$. To find $\epsilon y^{\prime}$, differentiate $y$ termwise, multiply by $\epsilon$ and then make a change of dummy variables $m=n-1.$ Then differentiate $y$ twice termwise and make a change of dummy variables $m=n-2$. The differential equation becomes $$\sum_{m=0}^{\infty}\epsilon^{m+2}a_{m+2}(m+2)(m+1)x^m+\epsilon^{m+2}a_{m+1}(m+1)x^m+\delta_{0,m}x^{m}=0$$ where $\delta_{0,m}$ is the kronecker delta. Putting $m=0,$we solve for $a_2.$ Then for $m>0$, we obtain a recursive relation. Then we can change back to the $y_n(x)$ notation.
A: 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.
