How to obtain $f(x)$, if it is known that $f(f(x))=x^2+x$? How to get $f(x)$, if we know that $f(f(x))=x^2+x$?
Is there an elementary function $f(x)$ that satisfies the equation?
 A: There is no such $f:\mathbb C\to\mathbb C$. See this paper:

When is $f(f(z)) = az^2 + bz + c$ ?
  by R. E. Rice, B. Schweizer and A. Sklar
The American Mathematical Monthly,
  vol. 87, no. 4 (Apr., 1980), pp. 252–263

More generally, they prove that a quadratic polynomial has no iterative roots of any order.
A: Firstly, let $g(x)$ be equal to $x^2+x$. Now we can say that $g(x)=g(y)$ is equivalent to $x=y$ or $x+y=-1$. So $g(g(x))=g(g(y))$ means that $g(x)=g(y)$ or $g(x)+g(y)=-1$. But $g(x)=\left(x+\frac{1}{2}\right)^2-\frac{1}{4}\ge-\frac{1}{4}$, so the second case can't take place. Thus, $g^n(x)=g^n(y)$ iff $g(x)=g(y)$ for every positive integer $n$.
Also $f$ returns each number not less than $-\frac{1}{4}$.
Well, I don't have a full solution and I print it on my little mobile phone, so my idea is the following: all real numbers can be divided onto many infinitive sequences, some of them are also infinitive to the left, every number is a value of $g$ at the previous number in the same sequence (if it exists), and some numbers don't occur in any sequence, but their values of $g$ are in one of our sequences. In other words, we draw an arrow from $x$ to $g(x)$ for every $x$. After all I've written and your imagination we can understand our situation. Now $f$ can just divide there sequences by pairs and map them to each other at every pair. Sorry for not too clear explanation
A: Well, one more answer, the most q&d trick, giving the formal power series for the half-iterate using the Newton-squareroot-algorithm applied to formal power series. (Thus it is in principle the same logic as the Carleman-ansatz as in my earlier answer  but looks stunningly simpler).
In Pari/GP one has the builtin-function "serreverse(f)" finding the reverse of a formal power series (which must not have a constant term as is in your problem).       
So we do the following     
   Z(x) = x + x^2       \\ define the function of which we want the half-iterate

        g = x + O(x^32)  \\ declare g as formal power series as initial "value"
   for(k = 2, 7, g = (Z(serreverse(g))+ g)/2 )    \\ just iterate several times

   print(g + O(x^9))      \\ correct to the seventh term:

Result: 
      x + 1/2*x^2 - 1/4*x^3 + 1/4*x^4 - 5/16*x^5 + 27/64*x^6 - 9/16*x^7 
       + 357/512*x^8 + O(x^9)

The eigth term were correct ( 171/256 x^8 )  if we had iterated one more time.

In the same way we would get the formal power series of the half-iterate for the notorious case $g(g(x)) = \exp(x)-1$ just by initializing $Z=exp(x)-1$ in the above code.
A: Here’s a technique for finding the first few terms of a formal power series representing the fractional iterate of a given function like $f(x)=x+x^2$. I repeat that this is a formal solution to the problem, and leaves unaddressed all considerations of convergence of the series answer.
I’m going to find the first six terms of $f^{\circ1/2}(x)$, the “half-th” iterate of $f$, out to the $x^5$-term. Let’s write down the iterates of $f$, starting with the zero-th.
\begin{align}
f^{\circ0}(x)&=x\\
f^{\circ1}=f&=x&+x^2\\
f^{\circ2}&=x&+2x^2&+2x^3&+x^4\\
f^{\circ3}&\equiv x&+3x^3&+6x^3& + 9x^4& + 10x^5& + 8x^6\\
f^{\circ4}&\equiv x &+ 4x^2& + 12x^3& + 30x^4& + 64x^5& + 118x^6\\
f^{\circ5}&\equiv x& + 5x^2& + 20x^3& + 70x^4& + 220x^5& + 630x^6\\
f^{\circ6}&\equiv x& + 6x^2& + 30x^3& + 135x^4& + 560x^5& + 2170x^6\\
f^{\circ7}&\equiv x& + 7x^2& + 42x^3& + 231x^4& + 1190x^5& + 5810x^6\,,
\end{align}
where the congruences are modulo all terms of degree $7$ and more.
Now look at the coefficients of the $x$-term: always $1$. Of the $x^2$-term? In $f^{\circ n}$, it’s $C_2(n)=n$. The coefficient of $x^3$ in $f^{\circ n}$ is $C_3(n)=n(n-1)=n^2-n$, as one can see by inspection. Now, a moment’s thought (well, maybe several moments’) tells you that $C_j(n)$, the coefficient of $x^j$ in $f^{\circ n}$, is a polynomial in $n$ of degree $j-1$. And a familiar technique of finite differences shows you that
\begin{align}
C_4(n)&=\frac{2n^3-5n^2+3n}2\\
C_5(n)&=\frac{3n^4-13n^3+18n^2-8n}3\,,
\end{align}
I won’t go into the details of that technique. The upshot is that, modulo terms of degree $6$ and higher, you have $f^{\circ n}(x)\equiv x+nx^2+(n^2-n)x^3+\frac12(2n^3-5n^2+3n)x^4+\frac13(3n^4-13n^3+18n^2-8n)x^5$.
Now, you just plug in $n=\frac12$ in this formula to get your desired series. And I’ll leave it to you to go one degree higher, using the iterates I’ve given you.
A: If we assume a function which has a power series expansion, we can just assume f to have a Maclaurin polynomial. Then the coefficients of $f(f(x))$ is a linear combination of iterated convolutions, for which the right hand side is $[0,1,1,0,0,\cdots]$. Also some regularization should likely be employed, punishing too large coefficients for large exponent monomials.
