The real part of $z=x+yi$ is given as $\Im f(z) = x^2-y_2+x$. Therefore I have $$ \left\{ \begin{array}{c} G_x = 2x+1 \\ G_y = -2y \end{array} \right. $$
$$ \left\{ \begin{array}{c} G_{xx} = 2 \\ G_{yy} = -2 \end{array} \right. $$ $G_{xx} + G_{yy} = 0$, therefore the function is harmonic.
Next step is $$ \left\{ \begin{array}{c} u_x = G_y = -2y \\ u_y = -G_x = -2x-1 \end{array} \right. $$
Then, $$u=\int u_x\ dx=-2xy+h(y)+C$$ $$u_y = -2x-1$$ Therefore, $$u = -2xy - 2x - 1$$ Finally, $$f(z) = u+iv = -2xy-2x-1 + i(x^2-y^2+x) +C= -2xy - 2x-1 + x^2i+y^2i+xi +C= i(x^2-y^2+2xyi) + (-2x-1+xi) + C = iz^2 + (-2x-1+xi) + C$$ But this is a wrong answer. Where is my mistake?
$u = -2xy - 2x - 1$ is wrong !
From $u=-2xy+h(y)+C$ we get
$u_y=-2x+h'(y)=-G_x=-2x-1$, hence $u=-2xy-y+c$.