Asymptotic development of a recurrent sequence Let $u_0 = 1$ and $u_{n+1} = \frac{u_n}{1+u_n^2}$ for all $n \in \mathbb{N}$.
I can show that $u_n \sim \frac{1}{\sqrt{2n}}$, but I would like one more term in the asymptotic development, something like $u_n = \frac{1}{\sqrt{2n}}+\frac{\alpha}{n\sqrt{n}} + o\bigl(\frac{1}{n^{3/2}}\bigr)$.
Here is the outline of my proof of $u_n \sim \frac{1}{\sqrt{2n}}$: 


*

*$(u_n)$ is decreasing, and bounded from below by $0$, hence converges.

*The limit $\ell$ satisifies $\ell = \frac{\ell}{1+\ell^2}$, hence $\ell = 0$.

*A computation gives $v_n = u_{n+1}^{-2} - u_n^{-2} \to 2$.

*Using Cesàro lemma, $\frac{1}{n} \sum_{k=0}^{n-1} v_k \to 2$.

*Hence $\frac{1}{n} u_n^{-2} \to 2$.

 A: Let $a_n=\dfrac1{u_n^2}$. Then
$$
a_{n+1}-a_n = \left(\frac{1+u_n^2}{u_n}\right)^2 - \frac1{u_n^2} = 2+u_n^2 = 2+\frac1{a_n}.
$$
From this and $a_2=4$ we can find a successive sequence of estimates for $n\ge2$:
\begin{align*}
 a_n &\ge 2n; \\
 a_n &\le 2n + \sum_{k=2}^{n-1}\frac1{2k}
      < 2n+\frac12\log n; \\
 a_n &\ge 2n + \sum_{k=2}^{n-1}\frac1{2k+\frac12\log k}
      > 2n+\int_2^n \frac{dx}{2x+\frac12\log x}
      \\ &> 2n+\int_2^n \frac{dx}{2x}
            - \int_2^n \frac{\frac12\log x}{2x(2x+\frac12\log x)}dx
      = 2n+\frac12\log n - \mathcal{O}(1).
\end{align*}
Then 
$$
u_n = a_n^{-1/2} =
\frac1{\sqrt{2n}} \left(1+\frac{\log n}{4n}+\mathcal{O}\bigg(\frac1n\bigg)\right)^{-1/2} \\
= \frac1{\sqrt{2n}} \left(1-\frac12 \cdot \frac{\log n}{4n} + \mathcal{O}\bigg(\frac1n\bigg) \right) =
\frac1{\sqrt{2n}} - \frac1{8\sqrt2}\cdot \frac{\log n}{n^{3/2}} + \mathcal{O}\left(\frac1{n^{3/2}}\right).
$$
A: Note: This answer does not contain a proof, but I'm fairly sure it's right. Posted in the hope that perhaps the answer will guide you towards a proof thereof.
Looking at the graph you can figure out that in fact the next correction looks to be $n^{-3/2} \log n$ dominating the next $n^{-3/2}$ term. Given this, one can do a plausible self-consistency check using the recurrence relation to find that the answer.
Specifically, suppose
$$u_n \sim \frac \alpha {\sqrt{2n}} + \frac {\beta \log n}{n^{3/2}}+\text{a sufficiently nice smaller series in particular }\mathcal{O}\left(\frac{1}{n^{3/2}}\right)$$
Then
$$u_{n+1}=\frac{u_n}{1+u_n^2}$$
implies
$$u_{n+1}-\frac{u_n}{1+u_n^2} \sim \frac 1 2 (\alpha - 2 \alpha^3) \frac 1 {n^{3/2}}+\cdots $$
so that $\alpha=0,-1/\sqrt 2,1/\sqrt 2$ are the options. You have checked that $\alpha=1/\sqrt{2}$.
Substituting this back in gives
$$u_{n+1}-\frac{u_n}{1+u_n^2} \sim \left(-\frac 1 {8\sqrt 2} - \beta\right) \frac 1 {n^{5/2}}+\cdots $$
and hence it seems the only consistent result has
$$\boxed{\displaystyle u_{n} \sim \frac{1}{\sqrt{2n}}-\frac{\log n}{8\sqrt 2 n^{3/2}}}$$
Numerically, I find the error
$$\epsilon = \frac{u_n - \frac{1}{\sqrt{2n}}}{-\frac{\log n}{8\sqrt 2 n^{3/2}}}  - 1\approx  \frac{3.4}{\log n} \to 0$$
suggesting that indeed this is correct with the next term in the series being the expected $n^{-3/2}$ term.
However, the coefficient of the $n^{-3/2}$ term cannot be determined by this self-consistency procedure, reflecting the fact that one can check this actually depends on the initial datum $u_0$.
Left to do: (Not that I think I'll come back to this myself, but others can!)


*

*Rigorously show the above logarithmic correction.

*Figure out the dependence of the next term on the initial data.



Mathematica code for my verification:
min = 2;
max = 5000000;
dat = RecurrenceTable[{u[n + 1] == u[n]/(N[1, 50] + u[n]^2), u[0] == 1}, u, {n, min, max}];
errordat = 
  Table[{M, (dat[[M - min + 1]] - 
        1/Sqrt[2 M])/(-Log[M]/(8 Sqrt[2] M^(3/2))) - 1}, {M, min, max, 100000}];
nextconst = ((dat[[max - min + 1]] - 
       1/Sqrt[2 max])/(-Log[max]/(8 Sqrt[2] max^(3/2))) - 1)*Log[max]
Show[ListPlot[errordat, AxesOrigin -> {0, 0}, PlotStyle -> PointSize[Large]], Plot[nextconst/Log[x], {x, min, max}, PlotStyle -> Red, PlotRange -> All]]

Output: Estimate of $3.446...$ for next constant, and a plot of the $\epsilon$ against $n$ with a fitted curve based on a logarithmic next correction:

Obviously the decay to 0 here is slow, which is to be expected given the predicted logarithmic correction. Feel free to go to larger ranges.
