<s>[INST] 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}
 [/INST] 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.
</s>