thoughts about $f(f(x))=e^x$

If $f$ is a solution, it is injective, hence increasing or decreasing. If $f$ was surjective, $\mathrm{exp}$ would be too, so $f$ is not. If $f$ was decreasing, it would have a fixed point, and $\mathrm{exp}$ too, so $f$ is increasing. The set $f(\mathbb{R})$ is equal to $]a,b[$ for some $a<b$, possibly infinite. If $b \neq + \infty$, $\exp (x) < b$ for all $x$, contradiction. So $b = + \infty$ and since $f$ is not surjective, $a \neq - \infty$. If $a \geq 0$, $\exp (x) \geq f(a)>0$ for all $x$, and that's impossible ($x=\log f(a)-1$ for example). So $a<0$. Since $\lim_{- \infty} f = a$, $\lim_{- \infty} \exp = f(a)$, and so $f(a)=0$.

Now the whole function $f$ can be reconstructed from $f|_{]- \infty, a]}$. Assume we only have a continuous, increasing and surjective function $f: ]-\infty,a] \rightarrow ]a,0]$. We can extend it to a continuous function on $\mathbb{R}$ such that $f \circ f = \exp$. For example, if $a < x \leq 0$, $x=f(y)$ for some unique $y \in ]-\infty,a]$, and so $f(x)=\exp y = \exp f^{-1}(x)$ ($f^{-1} : ]a,0] \rightarrow ]-\infty,a]$). It is easily checked that this defines a continuous, increasing and surjective function $f : ]-\infty,0] \rightarrow ]a,\exp a]$ (only the continuity at $a$ is not completely obvious, but easy). We can go on to extend $f$ on $]0,\exp (a)]$, $]\exp (a),1]$, etc.

So there are infinitely many continuous solutions, and a recipe that gives them all: pick $a<0$, and "pick" a continuous, increasing and surjective function $f : ]-\infty,a] \rightarrow ]a,0]$. There is a unique continuous extension of $f$ such that $f \circ f = \exp$.

If you want $f$ to be $C^k$, $C^{\infty}$, analytic; you need to take $f$ that is so on $]-\infty,a]$, and check that it is so around $a$ when you do the first extension (the behavior of $f$ on the right of $a$ is controlled by its behavior at $-\infty$, because $f(a+\epsilon)=\exp f^{-1} (a+\epsilon)$). Have fun!

This is a large and difficult topic. I put a number of references at LINK where one should begin with the obituary of Baker. Over the real line, the difficulty is cause by fixed points. So, while there are infinitely many complex numbers such that $e^z = z,$ for all real numbers $x$ we have $e^x > x$ and so $e^x \neq x.$ So my memory, and this should be in Baker's papers, and I believe is originally due to Kneser, is that there is a solution along the real line to $f(f(x)) = e^x,$ such that $f$ is real analytic, meaning it is analytic in a strip around the real line in the complex plane, but the "strip" may have varying width.

Sept. 2015: the actual reference is at volume 187 (1950):

Reelle analytische Lösungen der Gleichung $\varphi(\varphi(x)) = e^x$ und verwandter Funktionalgleichungen. Kneser, Hellmuth in: Journal für die reine und angewandte Mathematik | Journal für die reine und angewandte Mathematik | Article 56 - 67

Henning found a link

somehow I downloaded a copy from the Gottingen electronic library and have a pdf. This article is unusual, it concentrates on real analytic functions from the beginning (All of the Baker stuff is over $\mathbb C$), and there are no fixpoints of the function involved because we are restricting to the real line.

Note that the same question goes badly wrong if the target function decreases near a fixed point. So, there is a solution to $g(g(x)) = \sin x$ which is $C^\infty$ and actually analytic except possibly at integer multiples of $\pi/2.$ EDIT, September 2015: confirmed by J. Ecalle in email. In comparison, there is no such thing for $h(h(x)) = \cos x,$ as $\cos x$ is decreasing at the fixed point in $\mathbf R.$

The most flexible method for constructing these functions in the difficult case is apparently due to J. Ecalle. I do not have his papers. But the method is given in the book I call KCG,

Just one addendum to the present answers. Sometimes it seems to me, that the problem of non-uniqueness of the selection of some version of the half-iterate is taken too abstract. But in my view there should some natural argument, for instance some smoothness.
To show the effect of different initializations, which are so nicely worked out by @Plop , I made some pictures using excel. To have more information in the graphs I used however the function $ \small f(x)=h(h(x)) = b^x $ with $\small b=\exp(0.5)$ instead of $\small b=e $. These pictures are in (see: Excel_presentation) ; it is just a quick presentation; use only the first seven sheets. Here are three of that pictures so that you have an idea what I'm talking about:

It shows the functions $ \small f(x)=x, f(x)=b^x, f(x)=b^{b^x}, f(x)=b^{b^{b^x}} $ as strong curved lines, and the estimated 0,5,1.5,2.5 - iterated versions as far as computable from the initial guess (that a-value in Plop's post). For this the dottest lines were used to construct the subsequent values of the orbit. (If I recall right then this iterative construction was already used by G.H.Hardy, but I don't know whether he was the first one)
The goal of a good choice for a is then to make the magenta line for the half-iterate as smooth as possible.

But how to define/operationalize "smooth" with such a discrete orbit? One idea is to take the distances to the according coordinates of the mirrored function $\small f°^{-0.5}(x) $, see the next picture.

The distances show a sequence of erratically changing lengthes. This is shown in the next image.

We see, that a better choice of the initial value of a could provide a "smoother" version; instead of initializing with $\small (0,f° ^{0.5}(x))=(0,0.4) $ we use $\small (0,f° ^{0.5}(x))=(0,0.58) $ and get a far better curve of distances:

A set of 3 initial guesses is in the provided link of the Excel-html-output as given in the first link . The first three pictures are of the first type, then there is one picture explaining the distance-idea, and then the next three pictures are of the latter type. The following pictures are not of interest, they are just working material. You can use the tabstrip at the bottom of the image to navigate.