$f\colon\mathbb{R}\to\mathbb{R}$ such that $f(x)+f(f(x))=x^2$ for all $x$?

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.

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.