$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?

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!
