This (question 6) was from the November $2022$ NZ maths olympiad workshop (I couldn't find the online one, even on the Wayback Machine so I'm emailing them).
All I know is that $f$ is bijective (injective & surjective).
Assume that $f$ is not injective, ie there exist $x_1 < x_2$ such that $f(x_1) = f(x_2)$
$$
\begin{align}
f(x_1) &= f(x_2) \\[1ex]
f(f(x_1)) &= f(f(x_2)) \\[1ex]
x_1^3 + 1 &= x_2^3 + 1 \\[1ex]
x_1^3 &= x_2^3 \\[1ex]
x_1 &= x_2 \ \large{\unicode{x21af}}
\end{align}
$$
The last line is true since $x \mapsto x^3$ is an injective function. (it is strictly monotonic)
Assume that $f$ is not surjective, ie there exists a $y_0$ such that $\forall x : f(x) \neq y_0$
$$
\begin{align}
f(x) &\neq y_0 \\[1ex]
f(f(x)) &\neq f(y_0),\quad (f \text{ is injective}) \\[1ex]
x^3 + 1 &\neq f(y_0) \ \large{\unicode{x21af}}
\end{align}
$$
This is a contradiction since $x^3 + 1$ spans the reals, so $f$ is bijective. I have a feeling that $f$ is discontinuous, although I can't prove it.
I find it hard to even find specific values of $f$, like I only know that
$$f(\alpha) = \alpha$$
where $\alpha$ is the real solution to $x = x^3 + 1$. Trying to evaluate $f$ at other points results in other unknown values.
Below I have rearranged the functional equation (Wolfram Alpha couldn't solve it). Since $f$ is bijective there must be an inverse.
$$
\begin{align}
f(f(x)) &= x^3 + 1 \\[1ex]
f(x) &= f^{-1}\left(x^3 + 1 \right) \\[1ex]
f(x) &= \left[ f^{-1}(x) \right]^3 + 1 \ \text{(from eqn 1)}\\[1ex]
f^{-1}\left(x^3 + 1 \right) &= \left[ f^{-1}(x) \right]^3 + 1 \\[1ex]
\end{align}
$$
Alternatively one could rearrange the equation in terms of $f$
$$
\begin{align}
f(x) &= f^{-1}\left(x^3 + 1 \right) \\[1ex]
f\left( \sqrt[3]{x – 1} \right) &= f^{-1}(x) \\[1ex]
\left[ f\left( \sqrt[3]{x – 1} \right) \right]^3 + 1 &= \left[ f^{-1}(x) \right]^3 + 1 \\[1ex]
f(x) &= \left[ f\left( \sqrt[3]{x – 1} \right) \right]^3 + 1 \\[1ex]
f\left(x^3 + 1 \right) &= \left[ f(x) \right]^3 + 1
\end{align}
$$
Best Answer
This spells out an argument proposed by user tomasz.
First assume that $g:\mathbb R \to \mathbb R$ is 1-1, onto, and satisfies $g(x)>x$. Each orbit of the action of $g$ on $\mathbb R$ is a copy of $\mathbb Z$, with $g$ acting on each orbit by $n\mapsto n+1$. Let $X$ be the space of orbits. By choosing a single element in each orbit, we can consider that $X\subseteq \mathbb R$. Thus we have a 1-1 onto map $\phi\colon \mathbb R\to\mathbb Z\times X$, where $X\subseteq \mathbb R$. Via $\phi$, it suffices to define a suitable function $f$ ("operator-theoretic square root of $g$") on $\mathbb Z\times X$. Choose a fixed-point-free involution $\iota\colon X\to X$. We define $f$ by setting $f(n,x)=(n,\iota(x))$ if $x<\iota(x)$, and $f(n,x)=(n+1,\iota(x))$ otherwise. Then the composition $f\circ f$ sends $(n,x)$ to $(n+1,\iota\circ\iota (x))=(n+1,x)$, which corresponds to $g$ via the correspondence $\phi$.
Since each orbit is discrete, one can avoid the use of the axiom of choice in choosing a subset $X\subseteq \mathbb R$, by choosing the representative real number in each orbit to be the smallest positive number in the orbit. Doing the same for $\iota$ may require further work.
Now our $g(x)=x^3+1$ has a fixed point $-2<\alpha<-1$ so the proof will be slightly more complicated. First, we define $f(\alpha)=\alpha$. It remains to define $f$ separately on the intervals $(-\infty,\alpha)$ and $(\alpha,\infty)$. On the latter, $f$ satisfies the condition $f(x)>x$, whereas on the former, $f$ satisfies a similar condition $f(x)<x$, so a similar argument can be applied.