$f\colon\mathbb{R}\to\mathbb{R}$ such that $f(x)+f(f(x))=x^2$ for all $x$? A friend came up with this problem, and we and a few others tried to solve it.  It turned out to be really hard, so one of us asked his professor.  I came with him, and it took me, him and the professor about an hour to say anything interesting about it.
We figured out that for positive $x$, assuming $f$ exists and is differentiable, $f$ is monotonically increasing.  (Differentiating both sides gives $f'(x)*[\text{positive stuff}]=2x$).  So $f$ is invertible there.  We also figured out that f becomes arbitrarily large, and we guessed that it grows faster than any linear function.  Plugging in $f{-1}(x)$ for $x$ gives $x+f(x)=[f^{-1}(x)]^2$.  Since $f(x)$ grows faster than $x$, $f^{-1}$ grows slower and therefore $f(x)=[f^{-1}(x)]^2-x\le x^2$.
Unfortunately, that's about all we know...  No one knew how to deal with the $f(f(x))$ term.  We don't even know if the equation has a solution.  How can you solve this problem, and how do you deal with repeated function applications in general?
 A: Here is something: Iterating
$$f_0(x):=0, \quad f_{n+1}(x):= x^2-f_n\bigl(f_n(x)\bigr)$$
one obtains a sequence of polynomials whose lowest terms stabilize at
$$x^2-x^4+2x^6-4x^8+8x^{10}-16 x^{12}+32 x^{14}-65 x^{16}+\ldots-17316 x^{28} +\ldots\ .$$
The sequence of coefficients so obtained (apart from the sign) is listed here: OEIS, but I can't make anything out of the explanation given there.
A: I have a strategy to prove that there is such a function $f$ from $[0,2] \to [0,2]$. Intriguingly, it visually appears that this function may be unique!
Let $S$ be the square $[0,2] \times [0,2]$. Define $\sigma : S \to S$ by $\sigma: (x,y) \mapsto (\sqrt{x+y}, x)$. 
Let $f$ be any strictly increasing function $[0,2] \to [0,2]$ and let $\Gamma$ be the graph of $f$.
Lemma: The function $f$ obeys the required recursion if and only if $\sigma(\Gamma) = \Gamma$. 
Proof: Suppose that $f$ obeys the required recursion. We first check that $\sigma(\Gamma) \subseteq \Gamma$. Let $(x,f(x)) \in \Gamma$. Since $f$ is stricty increasing, there is some $x' \in [0,2]$ with $f(x') = x$. Then $f(x') + f(f(x'))= (x')^2$ or $x+f(x)=(x')^2$, so $x' = \sqrt{x+f(x)}$ and $f(x') = x$. We see that $\sigma(x,f(x)) = (x', f(x')) \in \Gamma$ as desired. 
We now check that $\Gamma \subseteq \sigma(\Gamma)$. Again, let $(x,f(x)) \in \Gamma$. Then we similarly have $\sigma(f(x), f(f(x)) )=(x,f(x))$.
Finally, suppose that $\sigma(\Gamma) = \Gamma$ for some increasing $f$. Then, for any $x \in [0,2]$, we would have $(x,f(x)) = \sigma(y, f(y))$ for some $y$. But then $y=f(x)$ and $x = \sqrt{y+f(y)}$, or $x = \sqrt{f(x)+f(f(x))}$ as desired. $\square$
So we need to find a path $\Gamma$, which is a graph of an increasing function, and is taken to itself by $\sigma$. A clearly necessary condition is that $\Gamma \subset \sigma^k(S)$ for every $k$.
Here is a plot of $\sigma^k(S)$ for $k=0$, $1$, ..., $6$.

Each $\sigma^k(S)$ is bounded between the graphs of two piecewise smooth functions, and they even appear to by Lipshcitz with a uniform constant for $k\geq 1$. My intuition is that, in that case, there always must exist a path between them. I haven't figured out how to make this precise yet; I'm posting in case someone else knows.
The really nifty thing is that it looks like $\bigcap_k \sigma^k(S)$ might itself be a path! In which case $f$ is unique! I 'm putting up this image now, to try to excite people into proving this.
A: A non-negative solution for $x \geq 0$ extends to all $x$ using $f(|x|)$.  Below I will consider only solutions of this type.
A non-negative continuous solution on $[0,+\infty)$ is injective and unbounded for positive $x$, and therefore increasing, with $f(\infty)=\infty$ a fixed point for $f$.  
Define an operator $G$ on functions by $G(f)=f(x)+f(f(x))$.
$f_0(x)=x^{\sqrt{2}}$ almost satisfies the equation, for large $x$, with $G(f_0) = x^2 + x^{\sqrt{2}} = x^2 + o(x^2)$.  
A correction $f_1(x) = x^\sqrt{2} - px^q$ can be found, with constants $p$ and $q$ chosen uniquely to cancel the $x^\sqrt{2}$ term so that $G(f_1)=x^2 + O(x^a)$ with $a < \sqrt{2}$, for large $x$.  The process can be repeated to get an asymptotic series with exponents in $\Bbb{Q}(\sqrt{2})$.  It is a reasonable guess that these are the asymptotics for large $x$, but proving the simplest case, that $\lim f(x)x^{-\sqrt{2}}=1$, or that there exists at least one solution with that behavior, is a challenging problem.  The simpler functional equation $F(F(x))=x^2$ has smooth solutions not asymptotic to $x^\sqrt{2}$.
If there are other fixed points of $f(x)$ they occur at $0$ or $2$.  
The series solution with fixed point at $x=0$ is $f(x)=x^2 - x^4 + 2x^6 - 4x^8 + 8x^{10} - 16x^{12} + 32x^{14} - \dots$ but unfortunately the next coefficients are $65, -138, 316$ and OEIS has it listed ( http://oeis.org/A141366 ) with data suggesting that the series does not increase very fast and might have nonzero radius of convergence (less than 1/(2.9)).  The ratio of adjacent coefficients increases from 2 to almost 3 over the first several hundred terms so of course it could also increase without bound.
The formal power series solution with fixed point $x=2$ has irrational coefficients.  It begins $f(2+t) = 2 + at + O(t^2)$ where $a= \frac{\sqrt{17}-1}{2}$ and the higher degree coefficients are in $\Bbb{Q}(a)$.
A: UPDATE: This answer now makes significantly weaker claims than it used to.
Define the sequence of functions $f_n$ recursively by
$$f_1(t)=t,\ f_2(t) = 3.8 + 1.75(t-3),\ f_k(t) = f_{k-2}(t)^2 - f_{k-1}(t)$$
The definition of $f_2$ is rigged so that $f_2(3) = 3.8$ and $f_2(3.8) = 3^2- 3.8 = 5.2$.
Set $g_k=f_k(3)=f_{k-1}(3.8)$. So the first few $g$'s are $3$, $3.8$, $5.2$, $9.24$, etc. Numerical data suggests that the $g$'s are increasing (checked for $k$ up to $40$). More specifically, it appears that 
$$g_n \approx e^{c 2^{n/2}} \quad \mbox{for}\ n \ \mbox {odd, where}\ c \approx 0.7397$$
$$g_n \approx e^{d 2^{n/2}} \quad \mbox{for}\ n \ \mbox {even, where}\ d \approx 0.6851$$  
We have constructed the $f$'s so that $f_k(3)=g_k$ and $f_k(3.8) = g_{k+1}$. Numerical data suggests also that $f_k$ is increasing on $[3,3.8]$ (checked for $k$ up to $20$). Assuming this is so, define
$$f(x) = f_{k+1} (f_{k}^{-1}(x)) \ \mathrm{where}\ x \in [g_k, g_{k+1}].$$
Note that we need the above numerical patterns to continue for this definition to make sense.
This gives a function with the desired properties on $[3,\infty)$.
Moreover, we can extend the definition downwards to $[2, \infty)$ by running the recursion backwards; setting $f_{-k}(t) = \sqrt{f_{-k+1}(t) + f_{-k+2}(t)}$. Note that, if $f_{-k+1}$ and $f_{-k+2}$ are increasing then this equation makes it clear that $f_{-k}$ is increasing. Also, if $g_{-k+1} < g_{-k+2}$ and $g_{-k+2} < g_{-k+3}$ then $g_{-k} = \sqrt{g_{-k+1} + g_{-k+2}} < \sqrt{g_{-k+2} + g_{-k+3}} = g_{-k+1}$, so the $g$'s remain monotone for $k$ negative. So this definition will extend our function successfully to the union of all the $[g_k, g_{k+1}]$'s, for $k$ negative and positive. This union is $[2, \infty)$.
There is nothing magic about the number $3.8$; numerical experimentation suggests that $g_1$ must be chosen in something like the interval $(3.6, 3.9)$ in order for the hypothesis to hold. 

I tried to make a similar argument to construct $f:[0,2] \to [0,2]$, finding some $u$ such that the recursively defined sequence $1$, $u$, $1^2-u$, $u^2-(1^2-u)$ etcetera would be decreasing. This rapidly exceeded my computational ability. I can say that, if there is such a $u$, then
$$0.66316953423239333 < u < 0.66316953423239335.$$
If you want to play with this, I would be delighted to hear of any results, but let me warn you to be very careful about round off errors!
A: Some observations (without assuming differentiability):
First, monotonicity in the left and right half line does not require differentiability.  Observe that if there exists $x_0, y_0\in\mathbb{R}$ such that $f(x_0) = f(y_0)$, then necessarily $f(f(x_0)) = f(f(y_0))$ and hence $x_0^2 = y_0^2$. So this implies that $f$ is injective among the positive (negative) numbers. If you assume we are looking for continuous solutions, this also implies that $f$ is monotonic in the region concerned.  Also, one gets that $f(x)$ cannot be bounded from above trivially because $x^2$ is not bounded from above. 
Next, observe that a fixed point of $f$ can only be $f(x) = x \implies 2x = x^2$, which means that the only possible fixed points of $f$ are $0$ and $2$. Along the same lines, we get that for a continuous solution, $f(0) \geq 0$: assume the contrary, then $f(0) = y < 0$. We have by the functional equation $f(y) = -y > 0$, so by continuity there exists some $y' < 0$ such that $f(y') = 0$. But this means that $f(y') + f(f(y')) = f(0) = y'^2 > 0$, a contradiction.   
Third, still assuming that $f$ is continuous, we ask whether $f$ can be bounded on either of the half lines. The answer is no, as it will necessarily contradict the functional equation. Hence $f(x)$ must be unbounded in each half line. Using that $f(x)$ must be unbounded from above, we have that there are three cases: $f(+\infty) < 0$, $f(-\infty) < 0$, or neither. (Cannot be both, because of monotonicity.) 
The second case can be ruled out, as for very large negative $M$, we would get $f(M) < 0$, so $f(f(M)) < f(0)$, and contradicting the functional equation. The third case will also be ruled out if $f(0) \neq 0$. By the previous arguments $f(0)$, if non-zero, must be positive, this implies that $f(f(0)) < 0$ and so $f$ cannot be monotonically increasing on the right line.  
In the first case, we get that monotonicity and unboundedness imply there exists $M > 0$ such that if $|x| > M$, $x f(x) < 0$ (in fact, using $f(0) + f(f(0)) = 0$, one can take $M= |f(0)|$). In particular this means that for $x < -M$, $f(f(x)) < f(0) < f(x) \implies f(x) \geq x^2 / 2$. But we can assume (by choosing a larger $M$ if necessary) that $M > 2$, which implies that $f(x) > M$, which implies that in fact $f(f(x)) < 0$ for $x < -M$. And hence $f(x) \geq x^2$ if $x < -M$. This implies that for sufficiently large and positive $x$, $x^2 = f(f(x))+ f(x) \geq f(x)^2 + |f(x)|$. This implies that for sufficiently large and positive $x$, $|f(x)| < x$.   
In the third case, monotonicity and positivity guarantees that $f(x) < x^2$ for $x > 0$. And given a solution on the right half line, setting $f(x) = f(-x)$ on the left half line gives automatically a continuous solution. Furthermore, we can show that $f(x)$ cannot be $O(x^\alpha)$ for any $\alpha < \sqrt{2}$. (Assume the contrary, for all sufficiently large $x$ we have $f(x) + f(f(x))\leq C x^{(\alpha^2)}$ for a universal constant $C$, and so contradicts the functional equation.) Similarly $f(x)$ cannot be bounded below, asymptotically, by any $\beta x^\alpha$ with $\alpha > \sqrt{2}$.  

BTW, I don't think your argument that "if $f$ is differentiable that $f(x)$ must be increasing for positive $x$" is correct. You used the fact that $f(x)$ is arbitrarily large (and positive) somewhere. But it doesn't have to be so for positive $x$: in that step you are making the assumption that $f:\mathbb{R}_+ \to \mathbb{R}_+$, which is not necessarily true. 
