Determine the general term of a Recurrent Rational Sequence How to determine the general term of the $(u_n)$ sequence with $u_0 \in \mathbb R_+$ and $u_{n+1}=\frac{u_n^2+1}{1+2u_n}$ ?
Source : les dattes à Dattier
 A: Part 1: An interesting phenomenon occurs for $u_0=2$:
$$u_1=1, u_2=\dfrac{2}{3}, u_3=\dfrac{13}{21}, u_4=\dfrac{610}{987}, u_5=\dfrac{1346269}{2178309}...$$
$$\text{i.e.,}
 \ \ u_1=\dfrac{F_1}{F_2}, u_2=\dfrac{F_3}{F_4}, u_3=\dfrac{F_7}{F_8}, u_4=\dfrac{F_{15}}{F_{16}}, u_5=\dfrac{F_{31}}{F_{32}}...$$
where the numerators and denominators are consecutive elements of Fibonacci sequence ($F_n$). See(https://oeis.org/A000045) with general form:
$$\tag{1}u_n=\dfrac{F_{2^n-1}}{F_{2^n}}$$ 
It is not so surprising because we know the connection between the Golden Ratio and Fibonacci sequence, with striking jumps (classical of the quick convergence of Newton's method ; see below). 
Relationship:
$$\dfrac{F_{2^{n+1}-1}}{F_{2^{n+1}}}=\frac{\left(\dfrac{F_{2^{n}-1}}{F_{2^{n}}}\right)^2+1}{1+2\dfrac{F_{2^{n}-1}}{F_{2^{n}}}}$$
is a consequence of identity $F_{2n-1}=F_{n}^2+F_{n-1}^2$ see this.

Part 2 (new Edit): Let us have a look at the general case now. We are going to transform the nice formula obtained by Francesco Alem., taking into account the fact that it has a narrow domain of validity (we must have $-1<\frac{2 u_0+1}{\sqrt{5}}<1$, which means $u_0<\dfrac{\sqrt{5}-1}{2}=\Phi-1.$)
Let us write this formula under the form:
$$\dfrac{2u_n+1}{\sqrt{5}}=\coth \left(2^n \coth ^{-1}\left(\frac{2 u_0+1}{\sqrt{5}}\right)\right)$$
or 
$$\tag{*}v_n=\coth \left(2^n \coth ^{-1}(v_0)\right)$$
by setting $$v_n:=\frac{2u_n+1}{\sqrt{5}}.$$
Let us transform $(*)$ by using the following formulas: 
$$\coth(x)=\dfrac{e^{2x}+1}{e^{2x}-1} \ \ \text{and} \ \ \coth^{-1}(x)=\frac12\ln\left(\dfrac{1+x}{1-x}\right)$$
(see (https://en.wikipedia.org/wiki/Inverse_hyperbolic_function)). Setting 
$$a=\dfrac{1+v_0}{1-v_0},$$
relationship (*) becomes
$$v_n=\left(exp(2  \ln(a^{2^{n-1}})+1)\right)/\left(exp(2 \ln(a^{2^{n-1}})-1)\right)$$

$$\tag{**}v_n=\frac{a^{2^{n}}+1}{a^{2^{n}}-1}$$

giving a very tractable explicit formula for $u_n$.
An important consequence of (**) is that we can in this way assert the convergence of sequence $v_n$, thus of sequence $u_n.$
The fact that $|a|>1$ implies that  $v_n\to1$. As a consequence, $u_n\to \frac{\sqrt{5}-1}{2}=\Phi-1.$
Remark : We have to look for another recurrence formula, complementary to (*) for the values of $v_0$ that are larger than $\Phi-1.$

Part 3: An interpretation of sequence $u_{n+1}=\frac{u_n^2+1}{1+2u_n}$.
It could come from Newton's method applied to the roots of $f(x)=x^2+x-1$:
$$u_{n+1}=u_n-\dfrac{f(u_n)}{f'(u_n)}=u_n-\dfrac{u_n^2+u_n-1}{2u_n+1}$$
that may converge towards one of the roots of $f$, i.e.,  $-\Phi$ or $\Phi-1$ (where $\Phi$ denotes the golden ratio: $\Phi\approx 1.618...$, i.e.) under the condition of being in one of the basins of attraction.
As $u_0>0$, convergence is necessarily towards the positive root $\Phi-1.$
"may converge" $\rightarrow$ means that a further analysis is necessary. The condition of convergence/divergence of Newton's method are not that evident...
A: 1/The trick :

Let $t,u,v$ functions with $t\circ v=v\circ u$, then : $t^n\circ v=v\circ u^n$.
Here, we have choosen, $t,u,v$ of the type : $$t(x)=\frac{x^2+1}{2x+1},
 u(x)=x^2 \text{ and } v(x)=\frac{ax+b}{cx+d}$$



2/For the solution proposed by Francesco, the form is :
$$h_n=f\circ g^n\circ f^{-1}$$

with $$f(x)=\frac{\sqrt 5}{2}\coth(x)-\frac{1}{2}$$
$$g(x)=2x$$
We have : $h_n\circ h_m=h_{n+m}$
and a conjecture: $$h_1(x)=\frac{1+x^2}{2x+1}$$

3/
we go to show : $h_1(x)=\frac{1+x^2}{2x+1}$
we have : $$(1):\coth(2x)=\frac{\coth(x)^2+1}{2\coth(x)}$$
More : $$f=w \circ \coth$$
with $$w(x)=\frac{\sqrt{5}}{2}x-\frac{1}{2}$$
and $$w^{-1}(x)=\frac{(2x+1)}{\sqrt{5}}$$
so 

$$h_1=w \circ \coth \circ g \circ \coth^{-1} \circ w^{-1} $$
$$\text{ with } w(x)=\frac{\sqrt{5}}{2}x-\frac{1}{2} \text{ and }  g(x)=2x$$

with using of $(1)$ then :
$$ \coth \circ g \circ \coth^{-1}(x)=\frac{x^2+1}{2x} $$
so 
$$h_1(x)=w\left(\frac{w^{-1}(x)^2+1}{2w^{-1}(x)}\right)=\frac{\sqrt{5}}{2}\left(\frac{(\frac{2x+1}{\sqrt{5}})^2+1}{2(\frac{2x+1}{\sqrt{5}})}\right)-\frac{1}{2}$$
$$h_1(x)=\frac{\sqrt{5}}{2}\left(\sqrt{5}\frac{(\frac{4x^2+4x+1}{5})+1}{2(2x+1)}\right)-\frac{1}{2}=\frac{\sqrt{5}}{2}\left(\frac{\sqrt{5}}{5}\frac{4x^2+4x+1+5}{2(2x+1)}\right)-\frac{1}{2}$$
$$h_1(x)=\frac{\sqrt{5}}{2}\left(\frac{1}{2\sqrt{5}}\frac{4x^2+4x+6}{(2x+1)}\right)-\frac{1}{2}$$
$$h_1(x)=\frac{1}{4}\left(\frac{4x^2+4x+6}{2x+1}-\frac{4x+2}{2x+1}\right)$$
$$h_1(x)=\frac{1}{4}\left(\frac{4x^2+4}{2x+1}\right)$$
$$h_1(x)=\frac{x^2+1}{2x+1}$$

4/ 
about the case $a_{n+1}=a_n^2-2$
This case is interested because, for the previous case, resolved by Francesco and Jean-Marie, we could make changes of invertible variables, here not.
We use, a function $v$, no invertible, so $u(x)=x^2$ and $t(x)=x^2-2$ are not conjugate, but the trick continues to work and here 
$$v(x)=x+\frac{1}{x}$$
we have : 

$$\text{ if } a_{n+1}=a_n^2-2, \text{ and }b \in \mathbb C \text{ with } a_0=b+\frac{1}{b} \text{ then } a_n=b^{2^{n}}+\frac{1}{b^{2^n}}$$

A: Mathematica has solved the problem.
$$
u_n=\frac{1}{2} \left(\sqrt{5} \coth \left(2^n \coth ^{-1}\left(\frac{2 u_0+1}{\sqrt{5}}\right)\right)-1\right)
$$
There is nothing intuitive about this solution imho.
Mathematica command:

FullSimplify[Part[RSolve[{a[n+1]==(a[n]^2+1)/(2a[n]+1),a[0]==s},a[n],n],2]]

