asymptotic expansion/approximation Find the small solution of $$y''-y\left ( 1-y^{2} \right )=0 \text{ with } y\left ( 0 \right )=\epsilon \ll 1$$
Making a pun, I decided that $$y^{3}\left ( 0 \right )\ll y\left ( 0 \right )$$ so neglect $$y^{3}\left ( 0 \right )$$
This gives $$y''_{0}-y_{0}=0$$
The general solution is $$c_{1}e^{-x}+c_{2}e^{x}$$ but I shall drop the non-decaying term to obtain $$c_{2}e^{-x}$$
$$y_{0}=\epsilon e^{-x}$$
obviously, the solution y has the form $$y=y_{0}+y_{1} $$where $$y_{0}\gg y_{1}$$
Substituting $$y=y_{0}+y_{1}$$ into the original equation I do not arrive at what my text states.

I desparately need a leg up to bridge this gap.
The whole numerical and asymptotic approach to me is extremely'fuzzy' 
Thanks in advance
 A: With this initial value and the "normal" look of the differential equation, one can expect that the solution staying small means that it is $O(ϵ)$ for all times $x\in[0,\infty)$.
Setting $y=ϵz$ one expects $z$ to stay bounded. $z$ now satisfies the scaled differential equation
$$
z''-z+ϵ^2z^3=0\text{ with } z(0)=1,
$$
which now has a clear perturbation term.
For the limit case $ϵ=0$, you get $z''-z=0$ with solution
$$
z(x)=ce^x+(1-c)e^{-x}
$$
For this to stay bounded we need $c=0$ so that the growing term vanishes.
For the next best approximation we look for a perturbation in the size of the perturbation parameter, $z=e^{-x}+ϵ^2z_1$. Insertion into the differential equation reduces to
$$
z_1''-z_1+(e^{-x}+ϵ^2z_1)^3=0.
$$
Again removing the small parts containing $ϵ$ gives 
$$
z_1''-z_1=-e^{-3x}\text{ with }z_1(0)=0.
$$
The bounded solutions of that have the form $z_1(x)=Ae^{-x}+Be^{-3x}$ which results in the equations for the coefficients
$$
(9B-B)e^{-3x}=-e^{-3x}\text{ and }A+B=0
$$
so that the new approximation is
$$
z(x)=(1+\tfrac18ϵ^2)e^{-x}-\tfrac18ϵ^2e^{-3x}+O(ϵ^4)
$$
or
$$
y(x)=ϵ(1+\tfrac18ϵ^2)e^{-x}-\tfrac18ϵ^3e^{-3x}+O(ϵ^5)
$$
A: We're going to substitute $y=y_0 + y_1$ into the differential equation then throw away most (but not too many) of the smallest terms. This last part is a bit tricky but it will start to make more sense as you become more familiar with it.
So, suppose that $y = y_0 + y_1$ where $y_1 \ll y_0$ (and $y_0$ is small as well). Substituting this into the equation
$$
0 = y'' - y(1-y^2)
$$
yields
$$
\begin{align}
0 &= (y_0 + y_1)'' - (y_0 + y_1) \Bigl[ 1 - (y_0 + y_1)^2 \Bigr] \\
&= (y_0 + y_1)'' - (y_0 + y_1) \Bigl[ 1 - y_0^2 - 2y_0y_1 - y_1^2 \Bigr]. \tag{1}
\end{align}
$$
Let's look inside of the brackets first, since there are a lot of very small terms there. The term $y_0^2$ is small since $y_0$ is small, but then $2y_0y_1$ is even smaller since $y_1 \ll y_0$, and thus $y_1^2$ is smaller still. We could throw all of these away, but then we would be left with the equation
$$
0 = (y_0 + y_1)'' - (y_0 + y_1)
$$
which is exactly the one we solved to find $y_0$. Since throwing all three of the small terms away doesn't give us any more information than we already have, let's try keeping the largest of them. Again, we know that
$$
y_0^2 \gg 2y_0y_1 \gg y_1^2,
$$
so let's only keep $y_0^2$. Equation $(1)$ is then approximated by
$$
0 \approx (y_0 + y_1)'' - (y_0 + y_1) \Bigl[ 1 - y_0^2 \Bigr], \tag{2}
$$
and now we simplify a little:
$$
\begin{align}
0 &\approx (y_0 + y_1)'' - (y_0 + y_1) \Bigl[ 1 - y_0^2 \Bigr] \\
&= y_0'' + y_1'' - y_0 - y_1 + y_0^3 + y_0^2y_1 \\
&= y_1'' - y_1 + y_0^3 + y_0^2y_1,
\end{align}
$$
where in the last line we used the fact that $y_0'' - y_0 = 0$.  Now $y_0^3$ is small but $y_0^2y_1$ is even smaller, so let's throw the smaller away. We then arrive at the final approximation to $(1)$,
$$
\begin{align}
0 &\approx y_1'' - y_1 + y_0^3 \\
&= y_1'' - y_1 + \epsilon^3 e^{-3x}. \tag{3}
\end{align}
$$
Solving this yields
$$
y_1(x) = -\frac{\epsilon^3}{8} e^{-3x} + Ae^{-x} + Be^x.
$$
We need to have $y_1 \ll y_0$, and this can't be true if $B \neq 0$, so we take $B = 0$. Further, if $A \neq 0$ then $y_1(x) \sim Ae^{-x}$ for large $x$, so we won't have $y_1 \ll y_0$ in that case either. Thus we take $A = B = 0$, and we are left with
$$
y_1(x) = -\frac{\epsilon^3}{8} e^{-3x},
$$
which is $\ll y_0$ both as $\epsilon \to 0$ and as $x \to \infty$.
Note that since $y = y_0 + y_1$ and $y(0) = y_0(0) = \epsilon$ we need to have $y_1(0) \ll \epsilon$, and this is indeed true since $y_1(0) = -\epsilon^3/8$.
Thus
$$
y(x) \approx \epsilon \left( e^{-x} - \frac{\epsilon^2}{8} e^{-3x} + \cdots \right).
$$
Note also that there is a typo in your exercise, the $+$ should be a $-$.
