Find $f(x)$ such that $f(f(x)) = x^2 - 2$

There is no such $f$.

From http://yaroslavvb.com/papers/rice-when.pdf , the question of existence is determined by the theorem:

Theorem 6. Let $\mathbb{R}$ be the real line. Let $g$ be a real quadratic polynomial, so that $$g(x)=ax^2+ (b + 1)x+c,$$ for all real $x$, where $a\ne 0$, $b$, and $c$ are in $\mathbb{R}$. ... set $\Delta(g)= b^2-4ac$. If $\Delta(g)> 1$, then g has no iterative roots of any order whatever. [That is, there is no $f$ such $f\circ f = g$.] If $\Delta(g) =1$, then $g$ can be embedded in a 2-sided flow on $\mathbb{R}$, all of whose members are continuous functions. If $\Delta(g) <1$, then $g$ can be embedded in a 1-sided flow on $\mathbb{R}$, all of whose members are continuous functions; but $g$ cannot be embedded in any 2-sided flow on $\mathbb{R}$.

As $\Delta(g) = 0 - 4(1)(-2) = 8 > 1$ in your case, the question of existence is negative.

Looking closely at the article, the main point is that no function with only one 2-cycle can have a square root. In our case that means that there can be no partial solution $f:D\to D$ of the funcional equation $f(f(x))=x^2-2$ in $D\subset\Bbb{R}$ if $x_0=\frac{-1+\sqrt{5}}{2}\approx 0.61803$ or $x_2=\frac{-1-\sqrt{5}}{2}\approx -1.61803$ are in $D$.

In fact, clearly $x_0^2-2=x_2$ and $x_2^2-2=x_0$ (this implies that $x_0\in D$ if and only if $x_2\in D$).

There can be no other pair $y_1\ne y_2$ with $y_1^2-2=y_2$ and $y_2^2-2=y_1$, since then $$ \{-1,2,x_0,x_2,y_1,y_2\} $$ would be roots of the polynomial $P(x)=x^4 - 4 x^2 - x + 2$, since $y_1^2-2=y_2$ and $y_2^2-2=y_1$ implies $$ (y_1^2-2)^2-2=y_1\quad\Rightarrow \quad y_1^4-4y_1^2+2=y_1 \quad\Rightarrow \quad P(y_1)=0 $$ and similarly $P(y_2)=0$.

Now, if $x_0\in D$ (or $x_2\in D$) and $f:D\to D$ satisfy $f(f(x))=x^2-2$, then $x_1:=f(x_0)$ and $x_3:=f(x_2)$ would be such a pair, a contradiction that proves $x_0\notin D$ (and $x_2\notin D$).

Hmm, if we look at $g(x) = f(f(x)) = x^2-2$ and find the fixpoints $t_0=-1,t_1=2$ then we might assume a solution $h(x)$ for $g(x) = h(x-2)+2 $ and then from this we search for the half-iterate of $h(x)$ defining $h(x-2)+2 = w(w(x-2))+2$ and shall have a solution of the problem by $f(x)=w(x+2)-2$, hopefully having nonzero range of convergence at least in the vicinity of the fixpoint $t_1$.

With this ansatz we find that $$h(x) = g(x+2)-2 = (x+2)^2-4 = 4x+1x^2$$ could be a replacement for $g(x)$ and from this we get the leading terms for the power series of its half-iterate $$w(x) = 2 x + 1/6 x^2 - 1/90 x^3 + 1/720 x^4 - 7/32400 x^5 \\ + 161/4276800 x^6 - 391/55598400 x^7 + O(x^8) $$ by standard algebraic manipulations. (To check this insert $w(x)$ for $x$ in that leading terms)

Thus $f(x)$ can be written as a power series around (fixpoint) $2$: $$f(x) = 2+ 2 (x-2) + 1/6 (x-2)^2 - 1/90 (x-2)^3 + 1/720 (x-2)^4 - 7/32400 (x-2)^5 \\ + 161/4276800 (x-2)^6 - 391/55598400 (x-2)^7 + O(x^8) $$ (The link gives a confirmation using WolframAlpha and a $\cosh()$ /$\text{arccosh}()$-formula, see my other answer)

By the first $64$ terms it looks as if $w(x)$ has a positive range of convergence; and computing a couple of examples gives $$ \begin{array} {c|c|c|c} x &g(x)=h(1-t_1)+t_1 & f(x) = w(x-t_1)+t_1 & f(f(x))=w(w(x-t_1))+t_1\\ \hline\\ 1 & -1 & 0.17942912356439677341 & -1.0 \\ 1.5 & 0.25 & 1.0431497617870845281 & 0.25 \\ 2.5 & 4.25 & 3.0403583699367069623& 4.25 \\ \end{array}$$ while the results for $f(x)$ are only (but well) approximated because I've only the truncation to say 32 or 64 terms so far and not yet an analytical description for the coefficients which only would allow an exact result.

More Explanations

1) Generalities (Remark: $f(x)$ is meant here in general, not yet your sought function)

There is consent (see for instance L. Comtet, "Advanced Combinatorics", pg 144-148) that for power series without constant term there is a meaningful fractional iterate using the Bell-polynomials based on the formal power series for the function. One can express this more elegantly and concise in terms of Bell/Carleman matrixes (see wikipeda but for convenience I use the Carleman matrix here in transposed form) , which have a form such that for a function $f(x)$ we have $$ V(x) \cdot B = V(f(x)) \\ V(x) \cdot B^2 = V(f(f(x))) \\ \cdots\\ V(x) \cdot B^n = V(f^{[n]}(x)) \\ $$ where the notation $V(x)=[1,x,x^2,x^3,...]$ and thus $V(f(x))=[1,f(x),f(x)^2,f(x)^3,...]$ means infinite vectors of an indeterminate argument $x$.

This notation allows to isolate the required coefficients of the power series of $f(x)$ into the (infinite, lower triangular) matrix $B$ . Then for the operation of composition and self-composition (aka iteration) via the Bell-polynomials it suffices to denote powers of the matrix $B$ - because they implicitely define that required Bell-polynomial driven operations on the formal powerseries .

If the power series for $f(x) = \sum_{k=1}^\infty b_k\cdot x^k$ has no constant term, then the associated Bell-/Carleman-matrix is lower triangular, and powers of it can be computed also on the finitely truncated versions - giving exact coefficients which are valid also for the untruncated/infinite-size case.
Moreover this is also possible for fractional powers, and in particular for the square-root, which provides then the solution that you are searching.

2) How this concerns your problem (here $f(x)$ means your function)

The procedure is now to find an equivalent power series (or polynomial) for your function (let's call it $g(x)=f(f(x))$) which has no constant term - by shifting the $x$ and $y$-coordinate, such that $g(x) = x^2-2 $ which has a constant term, is rewritten as $h(x) = 4x + x^2$ which has no constant term and then we have modeled for any number of iterations $$g(x) = h(x-2) + 2 \\ g(g(x)) = h(h(x-2))+2 \\ \cdots$$.

Then for the function $h(x)$ we build the (lower triangular) Carleman-matrix $H$ which begins as $$ H= \small \begin{bmatrix} 1 & . & . & . & . & . \\ 0 & 4 & . & . & . & . \\ 0 & 1 & 16 & . & . & . \\ 0 & 0 & 8 & 64 & . & . \\ 0 & 0 & 1 & 48 & 256 & . \\ 0 & 0 & 0 & 12 & 256 & 1024 \end{bmatrix} \qquad \text{ matrix H is of infinite size} $$ and obviously the dot product with a $V(x)$-vector evaluates then to a $V(h(x))$-vector: $$ V(x) \cdot H = V(h(x))$$

Then the squareroot (let's call it $W$) of $H$ such that $W^2 = H$ can be determined by solution of an infinite equation system or by diagonalization.
One solution is again lower triangular and instead of a polynomially function like $h(x)$ it provides the coefficients of a power series. It begins like $$ W = \small \begin{bmatrix} 1 & . & . & . & . & . \\ 0 & 2 & . & . & . & . \\ 0 & 1/6 & 4 & . & . & . \\ 0 & -1/90 & 2/3 & 8 & . & . \\ 0 & 1/720 & -1/60 & 2 & 16 & . \\ 0 & -7/32400 & 1/540 & 1/30 & 16/3 & 32 \end{bmatrix} \qquad \text{ matrix W is of infinite size}$$ (You can check it by hand when you compute the square of W)
The second column gives the coefficients for the power series of $w(x)$ with the property $w(w(x-2))+2 = g(x)$ and thus gives an acceptable solution $f(x)=w(x-2)+2$.

Additional remark: Often with fractional iterates the occuring power series have very small or zero range of convergence - then additional measures must be introduced to make the solution usable and meaningful, but the powerseries for $w(x)$ here seems to have a convergence-radius of about $4$ and can be used for a not-too-small range of $x$ -values.