Context:
I was working on answering this ODE question with a new method that I have never seen anyone use before (Frobenius method but with a twist) and I reduced the problem to a very difficult recurrence relation question; the terms in this case are the coefficients of a particular solution $y(x)$ with initial conditions $y(0)=0$ and $y'(0)=1.$
The recurrence relation is $c_{n\geq 3}=$
$$\frac1{n}\sum_{k=2}^{n-1}(n+1-k)c_{n+1-k}c_k-\frac{1}{n^2}$$
which I deduced is equivalent to
$$\frac{n+1}n\sum_{k=2}^{\lceil (n-1)/2\rceil}c_{n+1-k}c_k\\+\left(n\bmod 2\right)\frac{n+1}{2n} c_{(n+1)/2}^2\\-\frac1{n^2}$$
where $c_0=0, c_1=1, c_2=\frac14.$
The first few terms are
$$\left\{0,1,\frac14,-\frac5{72},-\frac{97}{1152},-\frac{10777}{172800},-\frac{243521}{6220800},-\frac{27545531}{1219276800},-\frac{1492739209}{78033715200},-\frac{212648420707}{21069103104000},-\frac{6885617878439}{1053455155200000},\cdots\right\}.$$
I was also able to deduce (though not rigorously so would be nice if someone could prove) that $\lim_{n\to\infty}c_n=0$ and $\sum_{n=0}^\infty c_n=1$ (which means the terms eventually become positive again and most probably alternate between the $2$ signs while becoming absolutely small to converge to $0.$)
I found (not rigorously) that for large $n\to\infty,$ the recurrence relation reduces to
$$c_{2n+1}=\pm c_{n+1}$$
which might be hinting at when the sign change happens.
I tried nonlinear substitutions like $c_n=e^n$ to extract the terms out of the sum, but I always get remaining terms in $n.$ I also tried $nc_n=d_n$ and then the same substitution $d_n=e^n$ but I get a constant equal to an $n$-dependent harmonic number. Other things I tried are expressing the terms in terms of $c_1$ (where $c_2=\frac{c_1}{2(3c_1-1)}$) for example to find a polynomial solution.
I tried the same stuff for the more general case of the ODE, where we have $c_0=y(0)≠0,$ but still couldn't solve it.
This is generally a problem when solving these types of non-linear ODEs, as you can see via my other answers, so it would be very helpful to have a general approach so that we can tackle the case when $c_\infty=\infty$ for example (where accurate calculation is necesssary,) though I think this is too ambitious.
I am not an expert in recurrence relations, so I would like the answers, if any, to be a bit comprehensive and beginner-friendly.
Any somewhat accurate approximating simpler sequence, insight on the sequence's behaviour; when the terms are positive/negative, are they bounded, (dis)proofs of what I deduced on my own heuristically (which might be very wrong,) how we could have maybe readjusted the original (general) problem to make it easier on ourselves to add more conditions in our toolkit like for example expanding $y(u=\frac1x)$ instead of $y(x)$ for $x\in(-\infty,\infty)\setminus [-1,1]\implies u\in(-1,1),$ which makes it easier to make an approximating series as the later coefficients don't matter as much, would be very helpful.
Thank you.
Best Answer
I have a problem but, for sure, I may be mistaken.
Starting from scratch, we need to solve $$y\,y''+y'(y'-1)-y=0 \qquad \text{with} \qquad y(0)=0 \,,\,\, \,y'(0)=1$$ with $$y=\sum_{n=1}^\infty c_n\,x^n$$ "My" first coefficients are $$\left\{1,\frac{1}{4},-\frac{1}{72},\frac{1}{288},-\frac{11}{10800 },\frac{337}{1036800},-\frac{16657}{152409600},\frac{92879}{24 38553600},-\frac{8971901}{658409472000},\frac{32752889}{658409 4720000},-\frac{17642030237}{9560105533440000},\frac{106092907 271}{152961688535040000},-\frac{330158692960783}{1252756229101 977600000},\frac{590867235322373}{5846195735809228800000},\cdots\right\}$$ which, for $2 \leq n \leq 200$ at least, are alternating in sign.
Using the coefficients for $10 \leq n \leq 200$, a linear regression gives, with $R^2>0.999999$, $$\log(|c_n|)=-0.2815-0.845614\, n-1.51134\, \log (n)$$
$$\begin{array}{l|lll} \text{} & \text{Estimate} & \text{Std Error} & \text{Confidence Interval} \\ \hline a & -0.2815005 & 0.0009018 & \{-0.2832795,-0.2797215\} \\ b & -0.8456144 & 3.7\times 10^{-6} & \{-0.8456216,-0.8456071\} \\ c & -1.5113390 & 0.0002813 & \{-1.5118940,-1.5107841\} \\ \end{array}$$
Concerning the absolute errors, the mean is $0.00065$ and the maximum is $0.00645$.
The sum of all coefficients (from $n=1$ to $n=200$) is $1.23880865563816$ which already atteined for $n=31$.
Edit
What I did by hand before this answer $$y=x+\sum_{n=2}^5 c_n\,x^n$$ leads to $$0=(4 a_2-1)x+(6 a_2^2-a_2+9 a_3)x^2+(20 a_2 a_3-a_3+16 a_4)x^3+(15 a_3^2+30 a_2 a_4-a_4+25 a_5)x^4+(42 a_3 a_4+42 a_2 a_5-a_5)x^5+ (28 a_4^2+56 a_3 a_5)x^6+72 a_4 a_5x^7+45 a_5^2x^8$$ giving $$a_2=\frac 14 \qquad a_3=-\frac{1}{72} \qquad a_4=\frac{1}{288} \qquad a_5=-\frac{11}{10800}$$ which are correct. What remains is $$-\frac{337 x^5}{28800}+\frac{1757 x^6}{1555200}-\frac{11 x^7}{43200}+\frac{121 x^8}{2592000}$$
Adding $a_6$ to the list gives $a_6=\frac{337}{1036800}$ (which is correct) and what remains is $$\frac{16657 x^6}{3110400}-\frac{601 x^7}{1036800}+\frac{4099 x^8}{27648000}-\frac{40777 x^9}{1119744000}+\frac{1249259 x^{10}}{179159040000}$$