Solution: To prove that \( u(x, y) \) is harmonic, we need to show that it satisfies Laplace's equation: \[ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0. \] At first, compute the partial derivatives of \( u(x, y) \) which gives \[ \frac{\partial u}{\partial x} = 3x^2 - 3y^2, \quad \frac{\partial u}{\partial y} = -6xy - 2. \] Now, compute the second partial derivatives, we get \[ \frac{\partial^2 u}{\partial x^2} = 6x, \quad \frac{\partial^2 u}{\partial y^2} = -6x. \] Adding these, \[ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 6x - 6x = 0. \] Thus, \( u(x, y) \) is harmonic.

Next, we find an analytic function \( f(z) = u(x, y) + \iota v(x, y) \), where \( z = x + \iota y \), such that \( \text{Re}(f) = u \). Recall that as $f$ is analytic, it must satisfy the Cauchy-Riemann equations, that is, \[ \frac{\partial u}{\partial x} = \frac{\partial v}{\partial y}, \quad \frac{\partial u}{\partial y} = -\frac{\partial v}{\partial x}. \] Using the Cauchy-Riemann equations: \[ \frac{\partial v}{\partial y} = 3x^2 - 3y^2, \quad \frac{\partial v}{\partial x} = 6xy + 2. \] Therefore, \begin{align*} v(x, y) = \int (3x^2 - 3y^2) \ \mathrm{d} y = 3x^2 y - y^3 + h(x), \end{align*} where \( h(x) \) is a function of \( x \) alone.

Differentiate \( v(x, y) \) with respect to \( x \) and equate to \( \frac{\partial v}{\partial x} \): \begin{align*} \frac{\partial v}{\partial x} = 6xy + h'(x) & \implies 6xy + 2 = 6xy + h'(x) \\ & \implies h'(x) = 2 \implies h(x) = 2x + C, \end{align*} where \( C \) is a constant. Therefore, \[ v(x, y) = 3x^2y - y^3 + 2x + C. \] So, the analytic function will be \begin{align*} f(z) & = u(x, y) + \iota v(x, y) \\ & = (x^3 - 3xy^2 - 2y) + \iota (3x^2y - y^3 + 2x + C). \end{align*} Using the condition \( f(0) = 0 \), we find \( C = 0 \). Therefore: \[ \textcolor{blue}{ \boxed{ f(z) = (x^3 - 3xy^2 - 2y) + \iota (3x^2y - y^3 + 2x) } } \]

Method 2: Milne-Thomson's Method

We could have used the Milne-Thomson's method to construct an analytic function \( f(z) \) from the real part \( u(x,y) \). Using Milne-Thomson's method, we construct the analytic function \( f(z) \) from the real part \( u(x, y) \) as follows: \begin{align*} f(z) & = u(z,0) - \iota \int \frac{\partial u}{\partial y} (z,0) \ \mathrm{d} z \\ & = \left( z^3 - 0 - 0 \right) - \iota \int (-0 - 2)\ \mathrm{d} z \\ & = z^3 + 2\iota z + C. \end{align*} Using the condition \( f(0) = 0 \), we get \( C = 0 \). Thus, the analytic function is \[ \textcolor{blue}{ \boxed{ f(z) = z^3 - 2\iota z } } \]