Q:

$f'=e^{f^{-1}}$, again

This question is a spin-off of this one, in which the OP asks whether there is a solution $f:\mathbb R\to\mathbb R$ of the functional equation (not exactly an ODE) $f'=e^{f^{-1}}$, where $f^{-1}$ is the compositional inverse of $f$. The posted answer exploits the growth of $f(x)$ when $x\to-\infty$ and obtains a contradiction, which resolves the question nicely, but also invites the following question: what if we restrict to $f:\mathbb R_{\ge0}\to\mathbb R_{\ge0}$ and impose $f(0)=0$? This idea has been explored in the comments, where a formal power series expansion is obtained for $f$ which does not seem to converge for any $x\ne0$.
Taking another approach, we can use an iteration scheme starting from $f_1(x)=x$ and inductively solve the ODE $f_{n+1}'=e^{f_n^{-1}}$ with the initial condition $f_{n+1}(0)=0$ to obtain $f_{n+1}$, much in the spirit of Picard iteration. Explicitly, for example, we have
$f_2'=e^x$ and $f_2=e^x-1$;
$f_3'=e^{\ln(x+1)}=1+x$ and $f_3=x+x^2/2$;
$f_4'=e^{\sqrt{1+2x}-1}$ and $f_4=e^{\sqrt{1+2x}-1}(\sqrt{1+2x}-1)$
and the next iteration produces non-elementary functions. It is clear that the sequence $(f_{2k-1})_{k\ge1}$ is increasing, $(f_{2k})_{k\ge1}$ is decreasing, and $f_{2k-1}<f_{2k}$, so there are respective limits $f_-=\lim_{k\to\infty} f_{2k-1}$ and $f_+=\lim_{k\to\infty} f_{2k+1}$, with $f_-\le f_+$. It is also clear that from $n\ge2$ on the function $f_n'=e^{f_{n-1}^{-1}}$ is positive and increasing, so $f_n$ is increasing and convex, which can be passed to the limit to show that both $f_-$ and $f_+$ are also increasing and convex. As such they are continuous, and by Dini's theorem $f_{2k-1}$ converges to $f_-$ locally uniformly and similarly for $f_+$. Furthermore, the inequality $|x-y|\le |f_n(x)-f_n(y)|$ (as $f_n'=e^{f_{n-1}^{-1}}\ge1$) can also be passed to the limit. Then the following chain of inequalities:
$|f_-^{-1}(x)-f_{2k-1}^{-1}(x)|\le |x-f_-(f_{2k-1}^{-1}(x))|=|f_{2k-1}(f_{2k-1}^{-1}(x))-f_-(f_{2k-1}^{-1}(x))|$
shows that $f_{2k-1}^{-1}$ converges locally uniformly to $f_-^{-1}$, which then implies $f_{2k}'$ converges locally uniformly to $e^{f_-^{-1}}$. Hence $f_+'=e^{f_-^{-1}}$, and similarly $f_-'=e^{f_-^{-1}}$. From this it can be shown that $f_{2k-1}$ converges to $f_-$ locally in $C^\infty$, so both $f_-$ and $f_+$ are smooth functions, and they form an orbit of order at most 2 of the above iteration scheme. Moreover it can be shown that the first $n$ terms of the Taylor expansion of $f_n$ agrees with what have been calculated formally in the previous comments, so both $f_-$ and $f_+$ have the same Taylor expansion as calculated using formal power series expansion.
In light of the above, a priori the following three scenarios can happen:

$f_-\neq f_+$ and we have a genuine orbit of order 2, consisting of two functions having the same Taylor expansion at 0 but not being identical.
$f_-=f_+$ is an actual solution to the equation $f'=e^{f^{-1}}$, but it is merely $C^\infty$ but not analytic, having a divergent power series expansion at 0.
$f_-=f_+$ is an actual solution to the equation $f'=e^{f^{-1}}$, and it is analytic on a neighborhood of 0; we are just misled by the first 100 or so terms of the Taylor expansion.

Now finally comes the question: which of the above scenario is the reality? In the first two scenarios, one can also ask what is the growth rate of $f_-(x)$ and $f_+(x)$ as $x\to+\infty$.

A:

There is no analytic local solution at $0$ to $f'=e^{f^{-1}}$, $f(0)=0$, that is, the formal power series solution is diverging. Together with the solution given in comments by fedja, this means the actual scenario is 2. For convenience of notation, I shall consider the equivalent equation 
$$\begin{cases} g' =e^{g\circ g},  \\ g(0)=0, 
 \end{cases}$$
satisfied by $g(x):=-f^{-1}(-x)$ (Indeed, by the rule of the derivative of an inverse, $(f^{-1})'(x)={1\over f'(f^{-1}(x))}=e^{-f^{-1}(f^{-1}(x))}$ so that $g'(x) =e^{g(g(x))}$; see also Tom Copeland's previous answer here.)
Indeed, assume by contradiction the formal power series solution $x+{1\over2}x^2+{1\over2}x^3+{2\over3}x^4+\&c.$ to the above equation has a positive radius of convergence. Then, it extends uniquely by analytic continuation to a maximally-defined analytic function, still denoted $g$ (that is, defined on the largest positive interval $[0,a)$, for some $0<a\le+\infty$). 
Note that the Taylor series of $g$ at $0$ has non-negative coefficients. This follows immediately by induction, equating the  coefficients of $g'$ and $e^{g\circ g}$; incidentally, this series is the EGF of the positive integer sequence OEIS A214645, as also remarked here. As a consequence (check the details below), $g$ is totally monotonic on $[0,a)$; in particular $g'(x)>g'(0)=1$ and $g(x)>x$ for all $0<x<a$, and $g$ is invertible. 
Then observe that $\log( g'( g^{-1}(x))$  is a well-defined analytic function on the interval $g[0,a)$, and coincides with $g$ locally at $0$. By the maximality of $[0,a)$ we have thus $g[0,a)\subset[0,a)$, but, due to the inequality $g(x)>x$ on $(0,a)$, this inclusion is only possible if $a=+\infty$, so that $g$ is unbounded. On the other hand, arguing as in Christian Remling's previous answer, since $e^{-g(g(t))}g'(t)=1$ and $g(t)\ge t$, we have for any $x\ge0$
$$x=\int_0^{x}e^{-g(g(t))}g'(t)dt=\int_0^{g(x)}e^{-g(s)}ds\le \int_0^{+\infty}e^{-s}ds=1 ,$$
a contradiction. 
$$*$$
Rmk 1. To justify the total monotonicity of $g$, note that, as a general elementary fact, a  real analytic function on an interval $I$, whose Taylor series at some point $x_0\in I$ has non-negative coefficients, has Taylor series with non-negative coefficients ay any point $x\in I$, $x\ge x_0$. Indeed, this is clear for $x_1\ge x_0$ within the radius of convergence of $x_0$, and since there is a uniform radius of convergence at any $y\in [x_0,x]$, one reaches $x$ by finitely many steps $x_0<x_1<\dots<x_n=x$. 
$$*$$
Rmk 2. The very same argument works for other differential-functional equations like e.g.
$$\begin{cases} g' =1 + {g\circ g},  \\ g(0)=0, 
 \end{cases}$$
that generates the sequence OEIS A001028. 
As before, a maximally-defined analytic solution $g$, if any, must be totally monotonic and defined for all  $x\ge0$, for otherwise $ g'\circ g^{-1} -1$ would be a proper extension of it. Then we reach a contradiction as before, with one more step needed: since we have ${  g'(t)\over 1+g(g(t))}=1$ and $g(t)\ge t$ for any $t\ge0$, we also have, for any $x\ge0$
$$x=\int_0^{x}{  g'(t)dt\over 1+g(g(t))}=\int_0^{g(x)}{  dt\over 1+g(t)}\le\int_0^{g(x)}{  dt\over 1+t}=\log(1+g(x)) ,$$
whence $e^x\le 1+ g(x)$; if we plug this into the latter inequalities again, we get
$$x=\int_0^{g(x)}{  dt\over 1+g(t)}\le \int_0^{g(x)}e^{-t}dt\le 1 ,$$
as before. By comparison, the same conclusion also holds for $g'=F( {g\circ g})$  with any $F$ analytic and totally monotonic on $(-\epsilon,+\infty)$, and with $F(0)=1$.

