Find three non-constant, pairwise unequal functions $f,g,h:\mathbb R\to \mathbb R$...

Here is a proof that no such continuous functions exist. I'll use juxtaposition to denote function composition. As noted in the comments, $$ h^2=fgh=f^2=ghf=g^2. $$ Let $e=f^4=g^4=h^4$. Note that $fgf=hf=g$ and $gfg=gh=f$, so $$ ef=g^4f=gf^2gf=gfg=f, $$ $$ fe=fg^4=fgf^2g=gfg=f. $$ Similarly $eg=g=ge$ and $eh=h=he$. In particular $e^2=ef^4=f^4=e$, so $e$ is idempotent. Consider its image $X=\mathrm{im}\,e\subseteq\mathbb R$, so $e|_X$ is the identity. We have $\mathrm{im}\,f\subseteq X$ since $ef=f$. Moreover $f|_X^{4}=e|_X=\mathrm{id}_X$, so $f|_X$ is a permutation of $X$, and similarly for $g$ and $h$. (aside: it now follows that $f|_X$, $g|_X$ and $h|_X$ generate a quotient of the quaternion group).

Suppose $|X|>1$. Since $e$ is continuous, $X$ is a (possibly infinite) interval. Any continuous permutation of such a set is strictly monotone (either increasing or decreasing). Moreover the only strictly increasing involution of $X$ is the identity. Indeed suppose $\sigma$ is such a function. If $\sigma(x)>x$ then $$ x=\sigma^2(x)>\sigma(x)>x, $$ a contradiction. Similarly we cannot have $\sigma(x)<x$, so $\sigma$ is the identity.

In particular $f|_X$, $g|_X$ and $h|_X$ are strictly monotone. Since $f|_X=g|_Xh|_X$, they can't all be decreasing. wlog suppose $f|_X$ is increasing. Then $f|_X^2$ is increasing. Since $f|_X^4=\mathrm{id}_X$, applying the above result twice gives $f|_X^2=\mathrm{id}_X$ and then $f|_X=\mathrm{id}_X$.

If $|X|=1$, then we also have $f|_X=\mathrm{id}_X$. In either case, for any $x\in\mathbb R$ we have $g(x)\in X$, so $$ h(x)=f(g(x))=g(x), $$ whence $h=g$.

(This is by no means a comprehensive and strict answer, just few observations I made).

If you removed $0$ from the domain, the functions might be:

$$f_0 (x) \equiv -x$$ $$g_0(x) \equiv \frac 1 x$$ $$h(x) \equiv -\frac 1 x$$

Or better don't remove $0$ but allow $\infty$ into the domain and the image, making each of them the projectively extended real line. This is interesting because if you take a broader look at complex numbers and their representation as the Riemann sphere, you will notice that:

• $f_0$ corresponds to rotating the sphere half turn, while the axis goes through $0$ and $\infty$;
• $g_0$ is a similar rotation, the axis goes through $-1$ and $1$;
• $h$ is also a similar rotation, the axis goes through $-i$ and $i$.

Rotating comes easy to my imagination so I will stick to this interpretation for a while, but we should remember that every half turn rotation is equivalent to some axial reflection.

So in this case the three functions (and their compositions) correspond to certain operations (and their compositions) in 3D space where we imagine the Riemann sphere is.

In general you could choose any three half turns with respect to mutually perpendicular axes. In this case however each of them must map the extended real line (which is a great circle on the Riemann sphere) onto itself. This means one of the axes must go through $-i$ and $i$, i.e. one of the functions must be our $h(x)$.

If I did my calculations right, more general forms are ($\alpha \in \mathbb R$):

$$f_\alpha (x) \equiv \frac {-x+\alpha} {\alpha x + 1},$$ $$g_\alpha (x) \equiv \frac {\alpha x + 1} {x-\alpha} $$ $$h(x) \equiv -\frac 1 x$$

The following picture shows the projectively extended real line as a cross section of the Riemann sphere. Few possible axes are drawn.

I said our functions correspond to certain reflections (or rotations) in 3D. Now in 2D the interpretation is:

• $f_\alpha$ and $g_\alpha$ correspond to reflections about certain axes;
• $h$ corresponds to the reflection through the point $S$ (or a half turn, if you wish).

To me the most surprising conclusion is that: $$f_0 (x) \equiv -x$$ and $$g_0(x) \equiv \frac 1 x$$ are more similar than I ever thought.

An example:

$$f(x) = \begin{cases} x, & \text{if $x\in \mathbb Q$} \\ -x, & \text{if $x\notin \mathbb Q$} \end{cases}$$

$$g(x) = \begin{cases} -x, & \text{if $x\in \mathbb Q$} \\ x, & \text{if $x\notin \mathbb Q$} \end{cases}$$

$$h(x)=-x$$

A simple answer that almost works is $f(x)=g(x)=-x,h(x)=x$ with the only problem being $f=g$ but we can patch that up by making each one the identity on different parts of the real line. Split the reals into $|x| \gt 2,1 \le |x| \le 2, |x| \lt 1$. Make each the identity on one part and $-x$ on the other two.

Explicitly $$f(x)=\begin {cases} x& x \lt -2\\-x & -2 \le x \le 2 \\x & x \gt 2 \end {cases}\\g(x)=\begin {cases} -x& x \lt -2\\x & -2 \le x \le -1 \\-x & -1 \lt x \lt 1\\x&1 \le x \le 2\\-x&2 \lt x \end {cases}\\h(x)=\begin {cases} -x& x \le -1\\x & -1 \lt x \lt 1 \\-x & x \gt 1 \end {cases}$$