If $f(x)=x^2-x-1$ and $f^n(x)=f(f(\cdots f(x)\cdots))$, find all $x$ for which $f^{3n}(x)$ converges. Let $f:\mathbb{R}\to\mathbb{R}$ be the polynomial defined by
$$f(x)=x^2-x-1$$
and let
$$g_0(x)=f(x),\quad g_1(x)=f(f(x)),\quad\ldots\quad g_n(x)=f(f(f(\cdots f(x)\cdots)))$$
The positive root of $f(x)$ is the famous golden ratio
$$\varphi=\frac{1+\sqrt{5}}{2}.$$

Problem: Find all $x\in\mathbb{R}$ for which the limit
  $$\lim_{n\to\infty}g_{3n}(x)$$
  exists.

The answer appear to be not at all trivial. Based on numerical computations, I conjecture that the limit exists if and only if
$$x\in[-\sqrt{2},1+\sqrt{2}]-\mathbb{N}.$$
But for some $x$, the convergence is very strange. For example, look at this plot of the sequence $\{g_{3n}(0.18)\}$ for $0\leq n\leq 200$:

EDIT: As pointed out by Zach Effman, there seem to be a fractal structure emerging from this problem. To better visualize it, I looked at the boundedness of the sequence $\{g_n(z)\}$ for $z\in\mathbb{C}$. More precisely, I made a grid of points $z\in\mathbb{C}$, colored black if $\{|g_n(z)|\}$ is bounded at $z$ and the more quickly $\{|g_n(z)|\}$ diverges to infinity the more white the pixel is. A dramatic fractal structure emerges:

Does anybody recognize this fractal? Is this a known fractal?
The center of this picture reveals a very nice structure:

 A: Here are some partial results. I show below that there are only three
possible limits (when the sequence is nonconstant), and for each limit I find an interval such that
if the initial value is in that interval then it converges to that limit.
Applying Zach Effman’s method to $g=f^3$ rather than to $f$, we have
$$
g(x)-x=(x^2-2x-1)(x^3-x^2-2x+1)^2 \tag{1}
$$
The roots of $x^2-2x-1$ are
$$
\alpha_1=1-\sqrt{2},\ \ \alpha_2=1+\sqrt{2} \tag{2}
$$
While the roots of $x^3-x^2-2x+1$ are
$$
\beta_1 \approx -1.246,\ \ \beta_2 \approx 0.445, \ \ \beta_3 \approx 1.801 \tag{3}
$$
We have
$$
g'(\alpha_1)=25-22\sqrt{2}\approx -6.11, \ \  
g'(\alpha_2)=25+22\sqrt{2}\approx 56.11 \tag{4}
$$
So both those fixed points are repulsive. On the other hand,
$$
g'(\beta_1)=g'(\beta_2)=g'(\beta_3)=1
$$
(because of the identity $g'(x)-1=(x^3-x^2-2x+1)(8x^4 - 20x^3 - 4x^2 + 18x + 2)$)
There is a unique $\gamma_1<\beta_1$ such that $g(\gamma_1)=\beta_1$. It is easy
to see then that the interval $I_1=[\gamma_1,\beta_1]$ is stable by $g$,
and if $x\in I_1$ then the sequence converges to $\beta_1$.
Similarly : there is a unique $\gamma_2\in (\beta_2,\beta_3)$ such that $g(\gamma_2)=\beta_2$, and then  the interval 
$I_2=[\beta_2,\gamma_2]$ is stable by $g$,
and if $x\in I_2$ then the sequence converges to $\beta_2$, and 
there is a $\gamma_3>\beta_3$ such that  the interval 
$I_3=[\beta_3,\gamma_3]$ is stable by $g$,
and if $x\in I_3$ then the sequence converges to $\beta_3$.
A: This answer is inspired by comments with Mercio.  First we map the problem to an equivalent $z \mapsto z^2+c\;$ Mandelbrot problem.  And then we can use all of the known tools for analyzing the Mandelbrot set to analyze the Op's problem.  For example, see Wolf Jung's mandelbrot link, which includes a Mandelbrot program with landing ray's.
$$z \mapsto z^2-\frac{7}{4}$$  
This is a conjugation of the original problem with z=x-1/2, and $g(x)=f(x)-1/2$  For clarity, call the new function $g(z)=z^2-7/4$.  We calculate the polynomial roots of  $g^{o3}(z)-z=0$, which factors into a product of two identical cubic equations, and a quadratic.  Here we show the parabolic 3-cyclic attracting fixed points, which are the roots of the cubic equation, followed by the roots of the 1-cyclic quadratic equation.
$$8z^3 + 4z^2 - 18z - 1=0\; \left(r_1= -1.7469796037; r_2=-0.054958132087; r_3= 1.3019377358\right)$$
$$4z^2 - 4z - 7\; z=\left(r_4= -0.914213562373, r_5=1.91421356237 \right)$$
If we iterate $z \mapsto g(z)$ an infinite number of times, every point that remains bounded is in the Julia set for $c=-7/4$.  Points that are completely inside a bug will eventually you get to a 3-cyclic attracting period 3-cycling towards the three parabolic fixed points mentioned above.  Otherwise, if there are point that escape to $\infty$ arbitrarily close to your point, then the point you picked is on the boundary of the Julia set.  And then we need to figure out the rotation angle for the landing ray which lands at your point, perhaps using the Bottcher function.
But the simplest way to handle points exactly on the Julia set boundary is to look at values of $g^{on}(z)$ for finite values of n.  If $g^{on}(z)=r_1,r_2,r_3, r_4, r_5$ for some finite value of n, then the 3-cycle converges, since once you land on one of these five points, iterating $g^{o3}(z)$ does not change the point. Other than these five cases, the function $g^{o3}(z)$ will not eventually repeat, nor will it attract towards a repeating point.  The three cycle parabolic attracting points are only attracting for points in a connected region; in other words inside a Julia set bug, and not on the boundary.
One can also use the Bottcher function to put iterating $z \mapsto z^2$ into correspondence with iterating $z \mapsto z^2-7/4$, as both functions iterate to infinity. There will be an analytic Bottcher function mapping $\infty$ to zero.  And then we generate the Bottcher function for $g(z)$
$$\beta\left(\frac{1}{g(1/z)}\right) = \beta(z)^2$$
The Bottcher function maps the unit circle to the boundary of the Julia set for $c=-7/4$.  So if we have a point on the boundary of the Julia set, then the point corresponds to a particular number in the Boettcher set of the form 
$$\exp^{2\pi i r}$$
Where r is the rotation angle of the point z.  We are particularly interested in these rotation angles:
$$r = \left( \frac{\pm 1}{7}, \frac{\pm 2}{7}, \frac{\pm 4}{7}, \frac{\pm 1}{3}, 1 \right)$$
Since landing rays with these three rotation angles correspond to the three parabolic fixed points mentioned above.  r=1 corresponds to the solution z=1.91421356237, and r=$\pm 1/3$ corresponds to -0.914213562373.  Notice that iterating $z \mapsto z^2$ will double the rotation angle r.  And if we start with r=1/7, then iterating $z^2-7/4$ will go through a 3-cycle iterating (1/7,2/7,4/7.  There are two landing rays for each of the parabolic cusps; so (-1/7,-2/7,-4/7) is identical.  In addition, we have the rotation angles for $r=(\pm 1/3, 1)$, which correspond to the two solutions for the primary fixed point.
If we look at the rational number for the rotation angle, $\frac{p}{q}$, for the point z, then we can tell how it will repeat.  In particular, if $q=7\cdot 2^n$ or $q=3\cdot2^n$, or $q=2^n$ then the point is a pre-periodic point, and after the appropriate number of iterations, $g^{o3}(n)$ will repeat exactly.  All other rotation angles will upon iterating $z \mapsto z^2 - 7/4$ will not go to one of these five points.
A: Not a complete answer, but this is what I can determine right away:
What you're describing is an iterated function. If the sequence $f^n(x)$ converges, then we must have $f(L) = L$ for the limit point. This gives $x^2-x-1=x$ with roots $1\pm\sqrt{2}$. At a fixed point, the absolute value of the derivative of our function, $|f'(x)|$ determines the local behavior of the iterated sequence. If the absolute value of the derivative is less than one, the fixed point is an attractor and starting points nearby converge to it. For the function here, $f'(x) = 2x-1$ giving $|f'(1-\sqrt{2})| = 2\sqrt{2}-1 \approx 1.8$ and $|f'(1+\sqrt{2}) = 1+2\sqrt{2} \approx 3.8$. Unfortunately, both are greater than one. This means the fixed points are repellors: starting points arbitrarily close to them are pushed away.
All I can say, then, is that the iterated function definitely converges for starting points $x_0 = 1\pm \sqrt{2}$ because in these cases the sequence is constant. Actually, that's not quite true: any sequence that hits those points will stay constant. So we could ask, what are the points $x$ such that $f(x) = 1 \pm \sqrt{2}$? They're $1\pm \sqrt{2}$ (of course) as well as $\pm\sqrt{2}$. Repeating this, these are the points I get:

I imagine there's some sort of fractal structure here.
