If $f^2$ is an analytic function then so is $f$

a) Since $f^2$ is holomorphic, the Cauchy-Riemann criterion (expressed in Wirtinger's wonderfully concise notation) $\partial f^2/\partial \bar z=2f\partial f/\partial \bar z=0$, shows that $f$ is holomorphic at all points $z$ where $f(z)\neq 0$, since there $\partial f/\partial \bar z=0$. [This is the part solved by Don Antonio in more classical notation]

b) Since $f^2$ is holomorphic its zeros are isolated and so are those of $f$ (they are the same!)
But $f$ is locally bounded at those potential singularities, because $f^2$ is, and so Riemann's theorem on removable singularities permits you to conclude that the singularity is bogus and that actually $f$ is holomorphic also at the zeros of $f$ .

Conclusion: $f$ is holomorphic everywhere on its domain.

Suppose

$$f(x,y)=u(x,y)+iv(x,y)=:u+iv\Longrightarrow f^2=u^2-v^2+2uvi$$

Since $\,f^2\,$ is analytic then the Cauchy-Riemann equations apply here:

$$(u^2-v^2)'_x=(2uv)'_y\;\;\;,\;\;\;(u^2-v^2)'_y=-(2uv)'_x$$

But

$$(u^2-v^2)'_x=2uu_x-2vv_x\;\;\,\;\;(2uv)'_y=2u_yv+2uv_y\\(u^2-v^2)'_y=2uu_y-2vv_y\;\;,\;\;-(2uv)'_x=-2u_xv-2uv_x$$

So equalling:

$$I\;\;\;\;uu_x-vv_x=\;\;\;\;uv_y+vu_y\Longleftrightarrow\;\, u(u_x-v_y)-v(u_y+v_x)=0\\II\;\;\;\;uu_y-vv_y\;=-uv_x-vu_x\Longleftrightarrow u(u_y+v_x)+v(u_x-v_y)=0$$

Can you take it now from here (i.e., check the CR equations for $\,f=u+iv\,$) ?