Asymptotic expansion of $u_{n + 1} = \frac12 \arctan(u_n)$ (I'm aware of Asymptotic expansion of $v_n = 2^nu_n$ where $u_{n+1} = \dfrac{1}{2}\arctan(u_n)$ but it has no answers…)
Let be $u_0 \in \mathbb{R}$ and the sequence $(u_n)_n$ defined by: $u_{n + 1} = \frac12 \arctan(u_n)$.
I define also: $v_n = 2^n u_n$, so I can show that: $\lim (u_n)_n = 0$ (by studying $x \mapsto \frac12 \arctan(x)$), thus, I can show that $(v_n)_n$ is monotone and converges because it is bound.
Now, I conclude: $u_n \sim \dfrac{l}{2^n}$, I'd like to determine $l$ more precisely.
Here is what I tried, I suspect $l$ to be something like $f(\pi)$ for some $f$ :


*

*push the asymptotic expansion of $\arctan$ to the 2nd order and reinject it ;

*use $\arctan(u_n) + \arctan(1/u_n) = \dfrac{\pi}{2}$ ;

*use series techniques to look for $\sum v_{n + 1} - v_n$, maybe conclude using Cesaro summation

 A: For convenience make a slight generalization of the problem. Let $\,u_0\,$ and $\,y\,$ be two given numbers and suppose
$\,u_{n+1} = y \arctan(u_n)\,$ for all $\,n\ge 0\,$ where $\ y=1/2\ $
in the original recursion. Define with power series the function
$$ F(x,y,z)\! :=\! z\!\left(\!x \!-\! \frac{1\!-\!z^2}{1\!-\!y^2}
\frac{x^3}3 \!+\!\frac{(1\!-\!z^2)((3\!-\!2z^2)\!+\!y^2(2\!-\!3z^2)}
{(1-y^2)(1-y^4)}\frac{x^5}{15}
   \!+\! O(x^7)\! \right) \tag{1}$$
which is required to satisfy the equation
$\,F(x,y,yz) = \arctan(F(x,y,z))\,y.\,$
Use iteration to get the equation
$\, u_n = F(x,y,y^n)\,$ where $\, x = \lim_{n\to\infty} u_n/y^n.\,$
More terms in the power series expansion can be easily found.
Equation $(1)$ implies
$\, u_n \approx y^n(x - (1-y^{2n})x^3/(3(1-y^2))).\,$
A: (For easier discussion, I suggest you to read the introduction of Schroder's equation and the section on 'Conjugacy' of iterated function, in case you are not familiar with these topics.)
Let $f(x)=\frac12\arctan x$, and $f_n(x)$ be the $n$th iteration of $f$.
Let us reduce functional iteration to multiplication: if we can solve the corresponding Schroder's equation
$$\Psi(f(x))=s\Psi(x)$$
then it is well known (and also straightforward) that 
$$f_n(x)=\Psi^{-1}(f'(a)^n\cdot\Psi(x))$$ where $a$ is a fixed point of $f$.

For the moment, let us focus on $\Psi(f(x))=s\Psi(x)$.
Clearly, in our case, $a=0$, and $s=f'(a)=\frac12$.

For $a = 0$, if $h$ is analytic on the unit disk, fixes $0$, and $0 < |h′(0)| < 1$, then Gabriel Koenigs showed in 1884 that there is an analytic (non-trivial) $\Psi$ satisfying Schröder's equation $\Psi(h(x))=s\Psi(x)$.

Thus, $\Psi$ is analytic.
A few more observations:


*

*$\Psi(0)=0$.

*$\Psi'(0)$ is up to our choice, since if a function $\psi$ is a solution to the Schroder's equation, then so is $k\cdot \psi$ for any constant $k$. For convenience, set $\Psi'(0)=1$.

*All other Taylor series coefficients of $\Psi$ are then uniquely determined, and can be found recursively. (The method will be illustrated below.)

*By Lagrange inversion theorem, $\Psi$ is invertible in a neighbourhood of $0$, and $\Psi^{-1}(z)=0+\frac1{\Psi'(0)}z+o(z)\implies \Psi^{-1}(z)\sim z\quad(z\to 0)$.

*Therefore, $f_n(x)=\Psi^{-1}(f'(a)^n\cdot\Psi(x))=\Psi^{-1}(2^{-n}\Psi(x))\sim 2^{-n}\Psi(x)$ as $n\to\infty$.


Hence, for the limit the OP wanted to evaluate, $$\ell:=\lim_{n\to\infty}2^nf_n(x_0)=\Psi(x_0)$$ 

We shall now determine all the Taylor series coefficients of $\Psi(x)$ (valid only for $|x|<1$), since it can be assumed $0\le x_0<1$. 
Obviously, $\Psi$ is an odd function. Let
$$\Psi(x)=x+\sum^\infty_{k=2}\phi_{2k-1} x^{2k-1}$$
The basic idea is to repeatedly differentiate both sides of $\Psi(f(x))=s\Psi(x)$ and substitute in $x=0$, then recursively solve for the coefficients.
For example, differentiating both sides three times and substitute in $x=0$, we obtain
$$-\Psi'(0)+\frac18\Psi'''(0)=\frac12\Psi'''(0)\implies\phi_3=-\frac49$$

Slightly modifying the notations of our respectable MSE user @Sangchul Lee, for $\lambda=(\lambda_1,\lambda_2,\cdots,\lambda_n)$ a $n$-tuple of non-negative integers:


*

*write $\lambda \vdash n$ if $\sum^n_{i=1}(2i-1)\lambda_i=2n-1$.

*write $|\lambda| = \sum_{i=1}^{n} \lambda_i$.

*define the tuple factorial as $\lambda !=\frac{|\lambda|!}{\lambda_1!\cdot\lambda_2!\cdots\lambda_n !}$.


I will state, without proof, Faà di Bruno's formula for odd inner function:
$$(\Psi\circ f)^{(2n-1)}=(2n-1)!\sum_{\lambda \vdash n}\lambda!\cdot\phi_{|\lambda|}\prod^n_{i=1}\left(\frac{f^{(2i-1)}(0)}{(2i-1)!}\right)^{\lambda_i}$$ 
$$\implies \frac12\phi_{2n-1}=\sum_{\lambda \vdash n}\lambda!\cdot\phi_{|\lambda|}\prod^n_{i=1}\left(\frac{(-1)^{i+1}}{2(2i-1)}\right)^{\lambda_i}$$ 

Further simplifications lead to the final result:
$$\ell=\Psi(x_0)=\sum^\infty_{k=1}\phi_{2k-1} x_0^{2k-1} \qquad{\text{where}}\qquad \phi_1=1$$
$$\phi_{2n-1}=\frac{(-1)^{n}}{2^{-1}-2^{1-2n}}\sum_{\substack{\lambda \vdash n \\ \lambda_1\ne 2n-1}}\phi_{|\lambda|}\frac{\lambda! (-1)^{(|\lambda|+1)/2}}{2^{|\lambda|}}\prod^n_{i=1}\frac1{(2i-1)^{\lambda_i}}$$
Yeah, I know it’s ugly. But that’s the best we can obtain.
If anyone have a nice math software, please help me calculate the first few Taylor coefficients.
A: The iteration has the form $$u_{n+1}=a_1u_n+a_3u_n^3+...$$ As usual in such situations (See the answer in Convergence of $\sqrt{n}x_{n}$ where $x_{n+1} = \sin(x_{n})$ with citation of de Bruijn: "Asymptotic Methods ..."), one can try a Bernoulli-like approach and examine the dynamics of $u_n^{-2}$. There one finds
$$
\frac1{u_{n+1}^2}=\frac4{u_n^2(1-\frac13u_n^2+\frac15u_n^4\mp...)^2}
=\frac4{u_n^2}+\frac83-\frac4{15}u_n^2+O(u_n^4)\tag1
$$
Thus for a first approximation use $$x_{n+1}=4x_n+\frac83\iff x_{n+1}+\frac89=4(x_n+\frac89)$$ so that $$u_n^{-2}\sim x_n=4^n(x_0+\frac89)-\frac89.\tag2$$
This gives as first approximation
$$
u_n\sim \frac{2^{-n}u_0}{\sqrt{1+\frac89u_0^2(1-4^{-n})}}.\tag3
$$

For the next term use $v_n=(u_n^{-2}+\frac89)^{-1}$ and express (1) in terms of $v_n$
$$
\frac1{v_{n+1}}=\frac4{v_n}-\frac4{15}\frac{v_n}{1-\frac89v_n}+O(v_n^2) \tag4
$$
so that 
$$
\frac1{v_{n+1}}-\frac{4}{15^2}v_{n+1}=\frac4{v_n}-\frac1{15} v_n - \frac{1}{15^2}v_n+O(v_n^2)=4\left(\frac1{v_{n}}-\frac{4}{15^2}v_{n}\right)+O(v_n^2) \tag5
$$
and consequently
$$
\frac1{v_{n}}-\frac{4}{15^2}v_{n}=4^n\left(\frac1{v_{0}}-\frac{4}{15^2}v_{0}+O(v_0^2)\right) \tag6
$$
As $\frac1v-\frac{4}{15^2}v=\frac1v(1-\frac4{15^2}v^2)$, leaving out the second term adds an error $O(v_n^2)$ which is a small fraction of $O(v_0^2)$. Thus
$$
\frac1{u_n^2}+\frac89=\frac1{v_n}=4^n\left(\frac1{v_{0}}-\frac{4}{15^2}v_{0}+O(v_0^2)\right)=4^n\left(\frac1{u_0^2}+\frac89-\frac{4}{15^2}\frac{u_0^2}{1+\frac89u_0^2}+O(u_0^4)\right)\tag7
$$
so that the improved approximation is
$$
u_n=\frac{2^{-n}u_0}{\sqrt{1+\frac89u_0^2(1-4^{-n})-\frac{4}{25}\frac{u_0^4}{9+8u_0^2}+O(u_0^6)}} \tag8
$$
A: Partial answer for $u_0>0$, then 
$$u_{n+1}-\frac{u_n}{2}=\frac{1}{2}(\arctan{u_n}-u_n)<0$$
because $f(x)=\arctan{x}-x<0$ for positive $x$, thus 
$$0<u_{n+1}<\frac{u_n}{2}<u_n \tag{1}$$
Using MVT, $\exists z\in(u_{n+1},u_n)$ s.t.
$$u_{n+1}-u_n=\frac{1}{2}\left(\arctan{u_n}-\arctan{u_{n-1}}\right)=
\frac{1}{2}\frac{u_{n}-u_{n-1}}{z^2+1}$$
or (because $\color{red}{u_n-u_{n-1}<0}$)
$$\frac{1}{2}\cdot \frac{u_{n}-u_{n-1}}{u_{n+1}^2+1}<
u_{n+1}-u_n<
\frac{1}{2}\cdot \frac{u_{n}-u_{n-1}}{u_{n}^2+1}$$
or
$$\frac{u_{0}-u_{1}}{2^n}\prod\limits_{k=1}^n\frac{1}{u_{k+1}^2+1}<
u_{n+1}-u_n<
\frac{u_{0}-u_{1}}{2^n}\prod\limits_{k=1}^n\frac{1}{u_{k}^2+1}$$
Considering $u_{n+1}-u_n \sim -\frac{l}{2^{n+1}}$ then
$$\frac{u_{0}-u_{1}}{2^n}\prod\limits_{k=1}^n\frac{1}{u_{k+1}^2+1}>
\frac{l}{2^{n+1}}>
\frac{u_{0}-u_{1}}{2^n}\prod\limits_{k=1}^n\frac{1}{u_{k}^2+1}$$
or
$$\frac{2(u_{0}-u_{1})}{\prod\limits_{k=1}^n\left(u_{k+1}^2+1\right)}=
\frac{2(u_{0}-u_{1})\left(u_{1}^2+1\right)}{\prod\limits_{k=1}^{n+1}\left(u_{k}^2+1\right)}>
l>
\frac{2(u_{0}-u_{1})}{\prod\limits_{k=1}^n\left(u_{k}^2+1\right)}$$
or 
$$L_1>l>L_2$$
where 
$$L_2=\frac{2(u_{0}-u_{1})}{\lim\limits_{n\to\infty}\prod\limits_{k=1}^n(u_{k+1}^2+1)} \text{ and } 
L_1=L_2\left(u_{1}^2+1\right) \tag{2}$$
So, it looks like Robert (see the comments) was right, it depends on $u_0$.

Note: the following limit exists 
$$\lim\limits_{n\to\infty}\prod\limits_{k=1}^n(u_{k+1}^2+1)$$
because
$$0<\sum\limits_{k=1}\ln(u_{k+1}^2+1)<\sum\limits_{k=1}u_{k+1}^2<\infty$$
by ratio test from $(1)$.

The following code is computing $(2)$ but with a $\frac{1}{u_0}$ factor. You will notice a certain stability for $\frac{L_1}{2^n u_0 \cdot u_n}$ and $\frac{L_2}{2^n u_0 \cdot u_n}$ for various $u_0$ 
from math import atan
from math import pow

N = 300
U_0 = 190.0

u = []

it = U_0
u.append(it)

for i in range(1, N):
    it = 0.5 * atan(it)
    u.append(it)

val = 1.0
for i in range(1, N):
    val *= (u[i] * u[i] + 1.0)


L2 = (2.0 * (u[0] - u[1]) / val) / u[0]
L1 = L2 * (u[1] * u[1] + 1.0)
MID = (L1 + L2) / 2.0

print "limit L1 =",L1
print "limit L2 =",L2
print "limit MID =",MID

for i in range(N-100, N):
    Lp1 = L1 / pow(2, i)
    Lp2 = L2 / pow(2, i)
    MIDp = MID / pow(2, i)

    r1 = Lp1 / u[i]
    r2 = Lp2 / u[i]
    rMID = MIDp / u[i]

    print Lp2," vs ",u[i]," vs ",Lp1," --- ",MIDp
    print r2," vs ",r1," --- ",rMID

Try it here.
