I asked a question at Math.SE last year and later offered a bounty for it, but it remains unsolved even in the simplest case. So I finally decided to repost this case here:
Is it possible to express the following indefinite integral in elementary functions? $${\large\int}\sqrt{x+\sqrt{x+\sqrt{x+1}}}\ dx$$
The answer is 'no'. Making the substitution $$ x = \frac{(t-1)(t-5)(t^2+2t+5)}{16t^2}, $$ one finds $$ {\textstyle\sqrt{x+\sqrt{x+\sqrt{x+1}}}\,\mathrm{d}x} = \frac{(t^2-2t+5)(t^2-5)\sqrt{t^4{-}2t^2{-}40t+25}\ \mathrm{d}t}{32t^4}. $$ Denote the right hand side of the above equation by $\beta$. Now, setting $$ Q(t) = \frac{(25-25t-14t^2-3t^3+t^4)}{96t^3}, $$ one finds that $$ \beta = \mathrm{d}\left(Q(t)\sqrt{t^4{-}2t^2{-}40t+25}\right) + \frac{(10+15t-2t^2+t^3)\,\mathrm{d}t}{8t \sqrt{t^4-2t^2-40t+25}}\, $$ The second term on the right hand side is in normal form for differentials on the elliptic curve $s^2 = t^4-2t^2-40t+25$ that are odd with respect to the involution $(t,s)\mapsto(t,-s)$. Since this term represents a differential that has two poles of order $2$ (over the two points where $t=\infty$) and two poles of order $1$ (over the two points where $t=0$), an application of Liouville's Theorem (on integration in elementary terms, with the differential field taken to be the field of meromorphic functions on the elliptic curve) shows that this term is not integrable in elementary terms. Hence $\beta$ is not, and hence the original integrand is not.
I'm adding a separate answer for the general question that the OP asked, which settles the question in the negative for all $n>2$ (and gives an alternate proof for $n=3$ to the one I gave above).
Recall that the OP defined a sequence of algebraic functions $f_n$ by the rule $f_0(x) = 1$, $f_1(x) = \sqrt{x+1}$, and $f_{n+1}(x) = \sqrt{x + f_n(x)}$ for all $n\ge 1$. It was observed that $f_n$ has an elementary antiderivative for $n=0$, $1$, and $2$, and the problem was to determine whether $f_n$ has an elementary antiderivative for some $n>2$.
I am going to show that there is no elementary antiderivative of $f_n$ when $n>2$.
Assume $n>2$ (NB: This is important, because the argument below will not work for $n\le2$; the reader may enjoy finding where it breaks down), and let $K_n = {\mathbb C}\bigl(x,f_n(x)\bigr)$ be the elementary differential field generated by $x$ and $f_n(x)$. Then $K_n$ is the field of meromorphic functions on the normalization $\hat C_n$ of the algebraic curve $C_n$ defined by the minimal degree $y$-monic polynomial $P_n(x,y)$ that satisfies $P_n\bigl(x,f_n(x)\bigr) \equiv 0$. This minimal degree is $2^n$; for example, $P_2(x,y) = (y^2-x)^2-x-1$ and $P_3(x,y) = \bigl((y^2-x)^2-x\bigr)^2-x-1$, etc.
Since $P_{n+1}(x,y) = (P_n(x,y)+1)^2-x-1$ for $n\ge 1$ with $P_1(x,y)=y^2-x-1$, one sees, by applying the Eisenstein Criterion to $P_n(x,y)$ regarded as an element of $D[y]$ with $D$ being the integral domain ${\mathbb C}[x]$, that $P_n(x,y)$ is irreducible for all $n\ge 1$. Hence, $\hat C_n$ is connected.
It will be important in what follows to observe that $K_n$ has an involution $\iota$ that fixes $x$ and sends $f_n(x)$ to $-f_n(x)$; this is because $P_n(x,y)$ is an even polynomial in $y$. The fixed field of $\iota$ is ${\mathbb C}\bigl(x,\,f_n(x)^2\bigr)$, and the $(-1)$-eigenspace of $\iota$ is ${\mathbb C}\bigl(x,\,f_n(x)^2\bigr)f_n(x) = K_{n-1}{\cdot}f_n(x)$.
Now, the curve $C_n\subset \mathbb{CP}^2$ has only one point on the line at infinity, namely $[1,0,0]$, but the normalization $\hat C_n$ has $2^{n-1}$ points lying over this point. They can be parametrized as follows: First, establish the convention that $\sqrt{u}$ means the unique analytic function on the complex $u$-plane minus its negative axis and $0$ that satisfies $\sqrt1 = 1$ and $\bigl(\sqrt{u}\bigr)^2 = u$. Let $\epsilon = (\epsilon_1,\ldots,\epsilon_{n-1})$ be any sequence with ${\epsilon_k}^2=1$ and consider the sequence of functions $g^\epsilon_k(t)$ defined by the criteria $g^\epsilon_1(t) = \sqrt{1+t^2}$ and $g^\epsilon_{k+1}(t) = \sqrt{1+\epsilon_{n-k}t g^\epsilon_k(t)}$ for $1\le k < n$. Choose, as one may, a $\delta_n>0$ sufficiently small so that, when $t$ is complex and satisfies $|t|<\delta_n$, all of the functions $g^\epsilon_k$ are analytic when $|t|<\delta_n$. In particular, one finds an expansion $$ g^\epsilon_n(t) = 1+\tfrac12\epsilon_1\,t + \tfrac18(2\epsilon_1\epsilon_2-1)t^2 + O(t^3). $$
Also, it is easy to verify that the disk in $\mathbb{CP}^2$ defined by $$ [x,y,1] = [1,\ t g^\epsilon_n(t),\ t^2]\qquad\text{for}\quad |t|<\delta_n $$ is a nonsingular parametrization of a branch of $C_n$ in a neighborhood of the point $[1,0,0]$. In the normalization $\hat C_n$, this is then a local parametrization of a neighborhood of a point $p_\epsilon\in \hat C_n$. Obviously, this describes $2^{n-1}$ distinct points on $\hat C_n$.
When $x$ and $f_n$ are regarded as meromorphic functions on $\hat C_n$, it follows that there is a unique local coordinate chart $t_\epsilon:D_\epsilon\to D(0,\delta_n)\subset \mathbb{C}$ of an open disk $D_\epsilon\subset \hat C_n$ about $p_\epsilon$ such that $t_\epsilon(p_\epsilon)=0$ and on which one has formulae $$ x = \frac1{{t_\epsilon}^2} \quad\text{and}\quad f_n(x) = \frac{g^\epsilon_n(t_\epsilon)}{t_\epsilon} = \frac{1+\tfrac12\epsilon_1\ t_\epsilon +\tfrac18(2\epsilon_1\epsilon_2-1)\ {t_\epsilon}^2} {t_\epsilon} + O({t_\epsilon}^2). $$ In particular, it follows that $f_n(x)$, as a meromorphic function on $\hat C_n$, has polar divisor equal to the sum of the $p_\epsilon$ and hence has degree $2^{n-1}$. Of course, this implies that the zero divisor of $f_n(x)$ on $\hat C_n$ must be of degree $2^{n-1}$ as well.
Note that the functions $g^\epsilon_k$ satisfy $g^{-\epsilon}_k(-t) = g^{\epsilon}_k(t)$, where $-\epsilon = (-\epsilon_1,\ldots,-\epsilon_{n-1})$. This implies that $\iota(p_\epsilon) = p_{-\epsilon}$ and that $t_\epsilon\circ\iota = -t_{-\epsilon}$.
Now, the $2^{n-1}$ zeroes of $f_n(x)$ on $\hat C_n$ are distinct, for they are the zeros of the polynomial $q_n(x) = P_n(x,0) = (q_{n-1}+1)^2-x-1$, and the discriminant of $q_n$, being the resultant of $q_n$ and $q_n'$, is clearly an odd integer, and hence is not zero. Thus, $C_n$ is a branched double cover of $C_{n-1}$, branched exactly where $f_{n}$ has its zeros. This induces a branched cover $\pi_n:\hat C_n\to \hat C_{n-1}$ that is exactly the quotient of $\hat C_n$ by the involution $\iota$ (whose fixed points are where $f_n$ has its zeros). Since one then has the Riemann-Hurwitz formula $$ \chi(\hat C_n) = 2\chi(\hat C_{n-1}) - B_n = 2\chi(\hat C_{n-1}) - 2^{n-1}, $$ and $\chi(\hat C_1) = \chi(\hat C_2) = 2$, induction gives $\chi(\hat C_n) = (3{-}n)2^{n-1}$, so the genus of $\hat C_n$ is $(n{-}3) 2^{n-2} + 1$. (This won't actually be needed below, but it is interesting.)
The only poles of $x$ and $f_n(x)$ on $\hat C_n$ are the points $p_\epsilon$,
and computation using the above expansions shows that,
in a neighborhood of $p_\epsilon$, one has an expansion of the form
$$
f_n(x)\,\mathrm{d} x
- \mathrm{d}\left(f_n(x)\bigl(\tfrac12\ x + \tfrac16\ f_n(x)^2\bigr) \right)
= \left(\frac{ (1-\epsilon_1\epsilon_2) }
{4{t_\epsilon}^2}
+ O({t_\epsilon}^{-1})\right)\ \mathrm{d} t_\epsilon\ .
$$
Thus, the meromorphic differential $\eta$ on $\hat C_n$
defined by the left hand side of this equation has, at worst, double poles
at the points $p_\epsilon$ and no other poles.
Now, by Liouville's Theorem, $f_n$ has an elementary antiderivative if and only if $f_n(x)\ \mathrm{d} x$ and, hence, the form $\eta$ are expressible as finite linear combinations of exact differentials and log-exact differentials.
Thus, $f_n(x)$ has an elementary antiderivative if and only if $\eta$ is expressible in the form
$$
\eta = \mathrm{d} h + \sum_{i=1}^m c_i\,\frac{\mathrm{d} g_i}{g_i}
$$
for some $h,g_1,\cdots g_m\in K_n$ and some constants $c_1,\ldots,c_m$. Suppose that these exist. Since $\eta$ has, at worst, double poles at the $p_\epsilon$ and no other poles, it follows that $h$ must have, at worst, simple poles at the points $p_\epsilon$ and no other poles; in fact, $h$ is uniquely determined up to an additive constant because its expansion at $p_\epsilon$ in terms of $t_\epsilon$ must be of the form
$$
h = \frac{\epsilon_1\epsilon_2-1}{4t_\epsilon} + O(1).
$$
Moreover, because $\eta$ is odd with respect to $\iota$, it follows that $h$ (after adding a suitable constant if necessary) must also be odd with respect to $\iota$. This implies, in particular, that $h$ vanishes at each of the zeros of $f_n$ (which, by the argument above, are simple zeros). This implies that $h = r\,f_n$ for some $r\in K_{n-1}$ that has no poles and satisfies $r(p_\epsilon) = (\epsilon_1\epsilon_2-1)/4$ for each $\epsilon$. However, since $r$ has no poles and $\hat C_n$ is connected, it follows that $r$ is constant. Thus, it cannot take the two distinct values $0$ and $-1/2$, as the equation $r(p_\epsilon) = (\epsilon_1\epsilon_2-1)/4$ implies.
Thus, the desired $h$ does not exist, and $f_n$ cannot be integrated in elementary terms for any $n>2$.
