I've been trying to find all functions $f(z):\Bbb{C}\longrightarrow\Bbb{C}$ that are analytic and satisfy
$$(\operatorname{Re}(f(z)))^2 = \operatorname{Im}(f(z))$$
After plugging in $f(z)=x+iy$, I got $f(z)=x+ix^2$.
These kind of functions never satisfy Cauchy-Riemann equations, since $\mathcal u_x=1 \neq 0=v_y$.
Does it follow from this argumentation that no such analytic $f(z)$ exists?
You might have rushed your analysis a little bit. The condition you have given is equivalent to saying that $$ f(x+ iy) = u(x,y) + i u(x,y)^2$$ Computing the Cauchy-Riemann equations we have $$ \frac{\partial u}{\partial x} = 2u\frac{\partial u}{\partial y}$$ and $$ \frac{\partial u}{\partial y} = -2u\frac{\partial u}{\partial x}$$ Combining these two we have $$\frac{\partial u}{\partial x} = - 4 u^2 \frac{\partial u}{\partial x} \iff \frac{\partial u}{\partial x} \left(1 + 4u^2\right) = 0$$ Since $1 + 4u^2 > 0$ for all $u \in \mathbb{R}$ we must have $\frac{\partial u}{\partial x} = 0$. In similar fashion we may also compute $$ \frac{\partial u}{\partial y} = - 4u^2 \frac{\partial u}{\partial y} \iff\frac{\partial u}{\partial y} \left(1 + 4u^2\right) = 0 $$ From which we may conclude that $\frac{\partial u}{\partial x} = 0$. So $u(x,y)$ must be a constant function. Since we did not specify any specific value of the constant function $u(x,y)$ then it follows that it works for all values. Thus functions with a constant real part which satisfy the hypothesis you have given are all solutions.
Another kind of solution: Note that ${\rm Im}(f(z))\geq 0$ for all $z$. Hence if we put $g(z)=\exp(if(z))$, we have $|g(z)|=\exp(-{\rm Im}(f(z))\leq 1$ for all $z$. By Liouville's theorem, $g$ is constant, and $f$ is also a constant, of the form $c+ic^2$, $c\in \mathbb{R}$.
If as usual we set
$z = x + iy \tag{1}$
and
$f(z) = u(x, y) + iv(x, y), \tag{2}$
then the condition
$(\operatorname{Re}(f(z)))^2 = \operatorname{Im}(f(z)) \tag{3}$
may be written
$v = u^2;\tag{4}$
thus $f(z)$ becomes
$f(z) = u(x, y) + iu^2(x, y). \tag{5}$
Since $f(z)$ is holomorphic, $u(x, y)$ and $v(x, y) = u^2(x, y)$ are harmonic; that is,
$\nabla^2u = 0 \tag{6}$
and
$\nabla^2 u^2 = 0. \tag{7}$
We compute $\nabla^2u^2$; we have
$\nabla^2 u^2 = \nabla \cdot \nabla u^2$ $= \nabla \cdot (2u \nabla u) = 2\nabla u \cdot \nabla u + 2u \nabla^2u$ $= 2 \vert \nabla u \vert^2, \tag{8}$
using (6); in deriving (8), we have also used the standard identity
$\nabla \cdot (gX) = \nabla g \cdot X + g \nabla \cdot X, \tag{9}$
holding for differeniable scalar functions $g$ and vector fields $X$; for more, see https://en.m.wikipedia.org/wiki/Vector_calculus_identities. It follows from (7) and (8) that
$\vert \nabla u \vert^2 = 0, \tag{10}$
or
$\vert \nabla u \vert = 0, \tag{11}$
or
$\nabla u = 0; \tag{12}$
we conclude that $u(x, y)$ must be a real constant $c$; it follows then that $v(x, y) = c^2$ and, finally $f(z)$ is a complex constant of the form
$f(z) = c + ic^2, \tag{13}$
where $c \in \Bbb R$.
Nota Bene: The technique used above, viz. showing that $\nabla u = 0$, may be applied to other, similar questions as well; see for example the
Related Problem:
