No es realmente una solución, pero aquí hay una perspectiva diferente que era demasiado larga para un comentario:
Tenemos $$ f(x) = \sqrt{x+f(x+1)} $$ o equivalente $$ f(x)^2 =x+f(x+1) $$ lo cual nos da una pseudo-recursión para $f'(x)$: $$ f'(x) = \frac{1+f'(x+1)}{2f(x)} $$ Simplemente formalmente, obtenemos:\begin{eqnarray} f'(x) = \frac{1}{2f(x)} + \frac1{4f(x)f(x+1)} +\frac1{8f(x)f(x+1)f(x+2)}+\cdots = \sum_{n=0}^\infty \frac1{2^{n+1}}\prod_{k=0}^n\frac1{f(x+k)} \end{eqnarray} ¿Converge esta serie? Parece cierto que debiera hacerlo, y de hecho, dado que $f(x)\ge\sqrt{x}$, obtenemos que el límite superior para cada término es $$ \frac1{2^{n+1}}\prod_{k=0}^n\frac1{\sqrt{x+k}} $$ lo cual es claramente sumable (y por ende también tenemos convergencia uniforme en intervalos cerrados cuando $x>0$, validando que esta es realmente la derivada de $f$, no simplemente de forma formal).
Además, para la $m$-ésima derivada cuando $m>1$, tenemos $$ \sum_{j=0}^m\binom{m}{j} f^{(j)}(x)f^{(m-j)}(x) = f^{(m)}(x+1) $$ lo cual puede resolverse para $f^{(m)}$ tomando $$ f^{(m)}(x) =\frac{\sum_{j=1}^{m-1}\binom{m}{j} f^{(j)}(x) f^{(m-j)}(x) + f^{(m)}(x+1)}{2f(x)} $$ y similarmente nos da (al menos ¡de forma formal!)$$ f^{(m)}(x) =\sum_{n=0}^\infty \frac{\sum_{j=1}^{m-1}\binom{m}{j} f^{(j)}(x+n) f^{(m-j)}(x+n)}{2^{n+1}\prod_{k=0}^n f(x+n)} $$ Debería ser posible mostrar de manera similar la convergencia uniforme de estas series.
Ahora, observa que la fórmula $$ f(x+1) = f(x)^2 - x $$ nos permite calcular explícitamente $f(x+n)$ en términos de $f(x)$ para todos los $n\in \mathbb{N}$. Digamos $f(x+n) = g_n(f(x),x)$. Entonces $f(x+n) = f(x+n-1)^2 - x-n = g_{n-1}(f(x),x)^2 - x - n$, entonces tenemos la siguiente recursión explícita para $g_n$:$$ g_n(y,x) = \begin{cases}y & n=0\\ g_{n-1}(y)^2 - x - n & n>0 \end{cases} $$ Por lo tanto, nuestra fórmula de derivada realmente se convierte en una ecuación diferencial:$$ f'(x) = \sum_{n=0}^\infty \frac1{2^{n+1}}\prod_{k=0}^n\frac1{g_k(f(x),x)} $$ Esto no es extremadamente útil, ¡pero es ligeramente interesante! También implicará que $f$ es la función diferenciable única que satisface $f(x)^2 = x+f(x+1)$ con su valor inicial $f(0)=\sqrt{1+\sqrt{2+\sqrt{3+\cdots}}}$.
Actualización: Una especie de solución. La raíz es $$ r = - \sqrt{1 + \sqrt{2 + \sqrt{3 + \sqrt{4 + \cdots} - \sqrt{ 1 + \cdots}} - \sqrt{1 + \sqrt{2 + \cdots} - \sqrt{1+\cdots}}} - \sqrt{1 + \sqrt{2 + \sqrt{3 + \cdots} - \sqrt{1 + \cdots}} - \sqrt{1 + \sqrt{2 + \cdots} - \sqrt{1+\cdots}}}} $$ En caso de no estar claro cuál es el patrón, cada vez que veas $n+\cdots$ léelo recursivamente como $$n+\cdots=n+\sqrt{n+1 + \cdots} - \sqrt{1+\cdots}$$
Para ver cómo esto es una solución, se puede ver a partir de esta lectura recursiva:\begin{eqnarray} \sqrt{1+r + \sqrt{2+r + \cdots}} = \sqrt{1 - \sqrt{1 + \cdots} + \sqrt{2 +\cdots}} = -r \end{eqnarray} y por ende $$ f(r) = \sqrt{r + \sqrt{1+r + \sqrt{2+r+\cdots}}} = \sqrt{r - r} = 0 $$ Una forma recursiva de encontrar esta raíz es observar, usando $f(x)^2 = x + f(x+1)$, que se tiene $f(r+1) + r = 0$ o equivalentemente, $r+1 = 1-f(r+1)$. Por lo tanto $r+1$ es el punto fijo único de $x\to 1 - f(x)$. Por lo tanto, formalmente se tiene$$ r+1 = 1-f(1-f(1-\cdots)) $$ De hecho, el punto fijo es atractor; uno puede mostrar inductivamente que $f_n'(1) < 1$ para todo $n$, y por ende $f'(1)<1$. Usando el hecho de que $r<-1$ y que $f'$ es decreciente, se tiene$$ f'(r+1) = \frac{1+f'(r+2)}{2f(r+1)} = \frac{1+f'(r+2)}{-2r} \le \frac{1+f'(1)}{-2r} < \frac{1}{-r} < 1 $$ Así, la derivada de $1-f(x)$ en $x=r+1$ es menor que $1$ en valor absoluto, por lo que iterar converge a su punto fijo. Definiendo una secuencia $c_n$ recursivamente por $$ c_{n} = 1 - f_n(c_{n-1}) $$ con $c_0 = 1$, se obtiene una secuencia que converge a $r+1$. Los primeros términos:\begin{eqnarray} c_0 &=& 1\\ c_1 &=& 0 \\ c_2 &=& 1-\sqrt[4]{1+\sqrt2}\\ c_3 &=& 1-\sqrt{1-\sqrt[4]{1+\sqrt2}+\sqrt{2-\sqrt[4]{1+\sqrt2}+\sqrt{3-\sqrt[4]{1+\sqrt2}+\sqrt{4-\sqrt[4]{1+\sqrt2}}}}} \end{eqnarray} claramente esto se vuelve bastante complicado bastante rápido. Programé esta recursión en Haskell y parece converger bastante rápido. Después de $c_{65}$ en mi computadora, muestra los valores oscilando entre $-0.21103728351247164$ y $-0.21103728351247142$, por lo que supongo que la raíz es $$ r\approx -1.211037283512471\dots $$ lo cual concuerda con los cálculos numéricos de otros usuarios.
