If we have the following iteration relation $$u_{n+1}=\epsilon+u_n+au_n^2$$ where $0<\epsilon<<1$ and $a>0$. The initial condition is $u_0=0$. What will be the leading order?
I have tried the following: rewrite to get $$u_{n+1}-u_n=\epsilon+au_n^2$$ Thus, we have $$u_n=u_n-u_0=(u_{n}-u_{n-1})+(u_{n-1}-u_{n-2})+\cdots+(u_1-u_0)$$ $$=n\epsilon+a(u_0^2+u_1^2+\cdots+u_{n-1}^2)$$
To the first order in $\epsilon$, it is simply $n\epsilon$. Thus we will have $n\sim\epsilon^{-1}$ if we want $u_n$ to be of order $O(1)$. However, if we go to second order, we have
$$u_1=\epsilon$$ $$u_2=2\epsilon+a\epsilon^2$$ $$u_3=3\epsilon+a(\epsilon^2+(2\epsilon+a\epsilon^2)^2)=3\epsilon+a(1^2\epsilon^2 +2^2\epsilon^2+O(\epsilon^3))$$
Inductively we see that the second order in $\epsilon$ is $$a(1^2+\cdots+(n-1)^2)\epsilon^2\sim an^3\epsilon^2$$ which is in some sense even faster than the first order term if we want $u_n\sim O(1)$ because now we only need to iterate $n\sim \epsilon^{-2/3}$ times.
I have also tried the third order, the result is $n^5\epsilon^3$. It seems to be a pattern of $n^{2k-1}\epsilon^k$ when we expand to the $k$-th order. Is it true? If so, is there an elegant way of proving it instead of painstakingly expansion?