The really clever trick to find $f(z)$ when we know that $f = u+iv$ and we have formulas for $u(x,y)$ and $v(x,y)$ is the following: Put $f = u + iv$, substitute $y=0$ and $x=z$ and your expession for $f(z)$ magically appears.
For your example:
$$ u + iv = a x^3 - 3dx^2 y - 3axy^2 + dy^3 + i(3ax^2y - 3dxy^2 + ay^3 + dx^3 + K)$$
Put $y=0$ and $x = z$ to obtain:
$$f(z) = az^3 + idz^3 + iK = (a+id)z^3 + iK.$$
This method works, but only if you already know that $u$ and $v$ are conjugate harmonic functions.
(An alternative, but much more cumbersome metod is to put $x = \frac12(z + \bar z)$ and $y = \frac{1}{2i}(z-\bar z)$, expand and simplify.)
Why does the trick work?
Let $g(z) = (a+id)z^3 + iK$. By construction $g(z) = f(z)$ whenever $z$ is real. Hence, by the identity theorem for holomorphic functions $g(z) = f(z)$ everywhere.
Variant 1: Harmonic functions have the mean value property,
$$f(z) = \frac{1}{2\pi} \int_0^{2\pi} f(z + re^{i\varphi})\,d\varphi$$
if $f$ is harmonic in $\Omega$ and $\overline{D_r(z)} \subset \Omega$.
$u$ is an entire harmonic function, hence
$$u(0) = \frac{1}{2\pi}\int_0^{2\pi} u(e^{i\varphi})\,d\varphi = \frac{1}{2\pi}\int_0^{2\pi} \cos(\varphi)e^{\cos\varphi}\cos (\sin\varphi) - \sin(\varphi)e^{\cos\varphi}\sin(\sin\varphi)\,d\varphi.$$
Whether we write $z$ and $re^{i\varphi}$ or $(x,y)$ and $(r\cos \varphi, r\sin\varphi)$ is completely immaterial. The complex notation is just more convenient sometimes.
Variant 2: Consider the analytic function $f = u+iv$.
Since $u$ and $v$ are both real, and $d\varphi$ is also real, we have
$$\begin{align}
\int_0^{2\pi} u(\cos\varphi,\sin\varphi)\,d\varphi &= \operatorname{Re}\left(\int_0^{2\pi} u(\cos\varphi,\sin\varphi)\,d\varphi + i\int_0^{2\pi} v(\cos\varphi,\sin\varphi)\,d\varphi\right)\\
&= \operatorname{Re} \int_0^{2\pi} f(e^{i\varphi})\,d\varphi.
\end{align}$$
Now,
The path integral over some closed curve is zero, over an analytic function.
is not correct as stated. On the one hand, the closed curve must not wind around any point in the complement of the function's domain - but since we have an entire function, that is vacuously satisfied here. More pertinent in the case at hand is that the integral theorem concerns only integrals with respect to $dz$ (it's a theorem about holomorphic differential forms), but here the integrand is $f(z)\,d\varphi$, not $f(z)\,dz$. Thus Cauchy's integral theorem does not apply.
However, for integrals over a circle, we have a simple correspondence between $dz$ and $d\varphi$. If we parametrise the circle as $\gamma(\varphi) = z_0 + r e^{i\varphi}$, then we have
$$dz = \gamma'(\varphi)\,d\varphi = ire^{i\varphi}\,d\varphi = i(z-z_0)\,d\varphi,$$
so we get
$$\int_0^{2\pi} f(e^{i\varphi})\,d\varphi = \int_{\lvert z\rvert = 1} f(z)\frac{dz}{iz},$$
and we see that that leads to Cauchy's integral formula,
$$\frac{1}{i} \int_{\lvert z\rvert = 1} \frac{f(z)}{z}\,dz = 2\pi\: f(0).$$
Best Answer
This is perhaps a more general approach.
Define $v(x,y)=\int_{0}^{y}u_{x}(x,t)\,dt+\phi(x)$ to be a guess for $v$. (We are guessing like this as $v$ must satisfy the CR equations).
Where $\phi$ is some arbitrary function of $x$.
then $v_{y}=u_{x}$ and $v_{x}=\int_{0}^{y}u_{xx}(x,t)\,dt+\phi'(x)$.
Now by Harmonicity of $u$ we get
$$v_{x}=-\int_{0}^{y}u_{yy}(x,t)\,dt+\phi'(x)=u_{y}(x,0)-u_{y}(x,y)+\phi'(x)$$
So if CR equations are to be satisfied then $u_{y}(x,0)+\phi'(x)$ should be $0$.
So you can solve the ODE in above and set $\phi(x)=-\int_{0}^{x}u_{y}(t,0)dt$
Thus $v(x,y)=\int_{0}^{y}u_{x}(x,t)\,dt-\int_{0}^{x}u_{y}(t,0)dt$