What I'm looking for is a function $\phi(x)$ such that $\phi(\phi(x))=f(x)$, where $f(x)=x^2-1$.
I am aware of this stack exchange post:
Square root of a function (in the sense of composition)
and furthermore, the post that was mentioned within that post:
https://mathoverflow.net/questions/17605/how-to-solve-ffx-cosx/44727#44727
However, using the resources in the latter link, I'm not sure how to start to come to a conclusion. In the top answer in the latter link, two criteria must be met (quoted straight from the answer):
1) The superfunction(flow) of f(x) grows not faster than an exponent
2) Runge phenomenon does not appear.
I have no idea how to test #1; that lies beyond the scope of my mathematical knowledge, and looking at the link for the Runge Phenomenon in #2 seems understandable enough, and testing with the analytic solution provided, it seems to lead to that phenomenon, but I have no idea how to come to that conclusion for sure. The answer also mentions that there are ways to combat Runge Phenomenon, but that also lies beyond the scope of my mathematical knowledge.
Best Answer
In General: Suppose we have $f(\alpha) = \alpha,$ while $f'(\alpha) = s$ with $0 < s < 1,$ also $f'(x) > 0$ and $f(x) \neq x$ for $x \neq \alpha,$ $C < x < D,$ for constants $C < \alpha < D.$ The first thing we do is move the fixpoint to the origin, with $$ g(x) = f(x + \alpha) - \alpha. $$ Then $g'(0) = s.$ As a result, pages 122-123 in KCG or pages 46-47 in Alexander or pages 77-79 in Milnor, we can solve Schroeder's equation $$ \psi(g(x)) = s \psi(x) $$ with $$ \psi'(0) = 1 $$
by $$ \psi(x) = \lim_{n \rightarrow \infty} \frac{g^n(x)}{s^n}, $$ where $g^1(x) = g(x)$ and then $$ g^n(x) = g \left( g^{n-1}(x) \right) $$ From above, $$ s \psi(x) = \psi(g(x)) $$ so $$ \psi^{-1} (s \psi(x)) = g(x ) = f(x+\alpha) - \alpha. $$ Substitute $x \mapsto x - \alpha$ to get $$ \psi^{-1} (s \psi(x - \alpha)) = f(x) - \alpha. $$ At this point, we can use $\sqrt s$ to define the half-iterate, $$ h(x) = \alpha + \psi^{-1} \left( \sqrt s \; \psi(x - \alpha) \right) $$ We will need the evident $$ \psi \left( h(x) - \alpha \right) = \sqrt s \; \psi(x - \alpha) $$ We may now simply calculate $$ h(h(x)) = \alpha + \psi^{-1} \left( \sqrt s \; \psi(h(x) - \alpha) \right) $$ $$ h(h(x)) = \alpha + \psi^{-1} \left( \sqrt s \; \psi( \alpha + \psi^{-1} \left( \sqrt s \; \psi(x - \alpha) \right) - \alpha) \right) $$ $$ h(h(x)) = \alpha + \psi^{-1} \left( \sqrt s \; \psi( \psi^{-1} \left( \sqrt s \; \psi(x - \alpha) \right) ) \right) $$ $$ h(h(x)) = \alpha + \psi^{-1} \left( \sqrt s \; \left( \sqrt s \; \psi(x - \alpha) \right) \right) $$ $$ h(h(x)) = \alpha + \psi^{-1} \left( s \; \psi(x - \alpha) \right) $$ $$ h(h(x)) = \alpha + f(x) - \alpha = f(x). $$
==========================================================================
This one has the advantage that the fixpoint at $\phi = \frac{1 + \sqrt 5}{2}$ has the function's derivative greater than one. This means that there is a real analytic half iterate for, at least, $x > 0.$ All you need to do is solve the Schroeder equation around $\phi.$ I will need to look up some things, but the references you want are Alexander and KCG.
Oh: it is not difficult to prove that there can be no answer for this that extends to the entire complex plane, or even extends to the entire real line. That is, however, not the end of the story.
Alright, pages 46 and 47 in Alexander, A History of Complex Dynamics. We need to use the inverse function to get derivative below $1,$ so we define $$ h(x) = \sqrt {1+x}. $$ Or pages 122-123 in Kuczma, Choczewski, and Ger, Iterative Functional Equations.
I need some more time to make an actual computer program. I have done the more difficult case with derivative $1,$ including program, http://mathoverflow.net/questions/45608/formal-power-series-convergence/46765#46765 and later http://math.stackexchange.com/questions/911818/how-to-obtain-fx-if-it-is-known-that-ffx-x2x/912324#912324
Anyway, as I said, your problem is conceptually far simpler.
It is necessary to move the fixpoint to the origin, so $$ g(x) = \sqrt {\phi^2 + x \;} - \phi. $$ We will then get the Schroeder function $$ \psi(x) = \lim_{n \rightarrow \infty} (2 \phi)^n g^n(x), $$ where $g^1 (x) = g(x)$ and then $$g^n(x) = g (g^{n-1}(x))$$ Here is an example with $x=1,$ showing $\psi(1) \approx 0.88523$
===========================================================
==========================================================