Low convergence of a perturbation solution to a non-linear ODE I have the following first-order non-linear ODE:
$$
\frac{1}{12}y' + \frac{1}{40} \epsilon \left( y' \right)^3 = -1  \tag{1}\label{1}
$$
where $y=y(x), \enspace x\in \left[0,1\right], \enspace \epsilon \ll 1$ is the perturbation parameter, and Eq. \eqref{1} is subject to,
$$
y(1) = 0. \tag{2}\label{2} 
$$
I sought a perturbation solution in powers of $\epsilon$ expanding $y$ as follows,
$$
y = y_0 + y_1 \epsilon + y_2 \epsilon^2 + \cdots. \tag{3}\label{3} 
$$
Substituting \eqref{3} in \eqref{1} and \eqref{2} produces the following problems:
$ O(1): $
$$
\frac{1}{12}y_0'  = -1 \\ y_0(1) = 0
$$
$ O(\epsilon): $
$$
\frac{1}{12}y_1' + \frac{1}{40} \left( y_0' \right)^3 = 0 \\ y_1(1) = 0
$$
$ O(\epsilon^2): $
$$
\frac{1}{12}y_2' + \frac{3}{40} \left( y_0' \right)^2 y_1' = 0 \\ y_2(1) = 0
$$
The corresponding solutions are:
$$
y_0 = 12(1-x), \tag{4}\label{4} 
$$
$$
y_1 = -\frac{2592}{5}(1-x), \tag{5}\label{5} 
$$
$$
y_2 = \frac{1679616}{25}(1-x). \tag{6}\label{6} 
$$
Then, inserting Eqs. \eqref{4}-\eqref{6} in \eqref{3} yields the perturbation solution up to $ O(\epsilon^2) $. Finally, the exact solution of \eqref{1} is presented to measure the accuracy of \eqref{3}, namely
$$
y = \frac{ \sqrt[3]{100} \epsilon - \sqrt[3]{10 \alpha^2} }{ 3 \epsilon \sqrt[3]{\alpha} } (1-x), \tag{7}\label{7} 
$$
where $ \alpha = \sqrt{ 2 \epsilon^3 \left( 1458 \epsilon + 5 \right)  } - 54 \epsilon^2 $.
Now, the exact solution \eqref{7} is plotted with its corresponding numerical solution (via MATLAB ode15i function). Figure 1 shows a perfect agreement between the exact and the numerical solution.
Figure 1. Exact and numerical solutions
Here comes my problem (after all this wordiness). Figure 2 shows that the perturbation solution \eqref{3} only converges to the exact solution \eqref{7} for really small values of $\epsilon$, namely, $\epsilon \sim O(10^{-3})$. The title in each subplot shows the $\epsilon$ value and the relative error of the perturbation solution to the exact solution. This means that the contributions of $O(\epsilon)$ and higher-order terms are not important, and all the physics of the problem is driven by the leading-order solution $y_0$ \eqref{4}.
Figure 2. Perturbation and exact solutions
Questions:

*

*Did I achieve a correct perturbation solution?

*Is this non-linear equation not suitable for a regular expansion approach?

*Why does the perturbation parameter need to be so small?

A final comment. I proposed $ y = y_0 + y_1 \epsilon^\beta + \cdots,$ and I got $ \beta = 1$.
 A: Writing the solution as
$$
y(x) = \gamma(\epsilon) (1-x)  \ \ \ \quad \mathrm(0)
$$
we see that $\gamma$ must satisfy
$$
\gamma^3\frac{\epsilon}{40}+\frac{\gamma}{12} -1=0.  \ \ \ \quad \mathrm(1)
$$
This is a cubic equation so it always has one real solution which is
$$
g(\epsilon) = \frac{1}{3} \left(\frac{10^{2/3}}{\sqrt[3]{\sqrt{2}
   \sqrt{1458 \epsilon ^4+5 \epsilon ^3}-54 \epsilon
   ^2}}-\frac{\sqrt[3]{10} \sqrt[3]{\sqrt{2}
   \sqrt{1458 \epsilon ^4+5 \epsilon ^3}-54 \epsilon
   ^2}}{\epsilon }\right).
$$
Now, the polynomial in (1) changes order when $\epsilon=0$ so you don't even have guarantee that the solution is continuous in $\epsilon$ at zero. Indeed $g(\epsilon)$ is complex for $\epsilon<0$. Computing the derivatives as $\epsilon \to 0^+$ one obtains
$$
g(\epsilon) = 12-\frac{2592 \epsilon }{5}+\frac{1679616 \epsilon
   ^2}{25}+O\left(\epsilon ^{3}\right),  \ \ \ \quad \mathrm(3)
$$
so your calculation is correct. However  the coefficients $a_n$ of this expansion defined as
$$
g(\epsilon) = \sum_{n=0} a_n \epsilon^n   \ \ \ \quad \mathrm(4)
$$
grow exponentially $a_n \sim b a^n$ (with $a\approx 186.9$, $b\approx3.14$) which means that the series (4) diverges for $\epsilon a >1$.  I think this explains what you are observing. You could make it even more evident by plotting the exact form of $g(\epsilon)$ together with one of its Taylor approximants.
