Finding general harmonic polynomial of form $ax^3+bx^2y+cxy^2+dy^3$.

An alternative approach is to 'guess' the form of the final analytic function. We want $f(z)=u(x,y)+i v(x,y)$ where $z=x+i y$ and $u(x,y)$ is a homogeneous cubic polynomial. Since $u(x,y)$ is cubic, we infer that $f(z)$ should also be some cubic polynomial; since it is homogeneous, we further conclude that $f(z)=\alpha z^3$ for some $\alpha$. (If it contained $z^2$, for instance, $u(x,y)$ would have terms of degree 2 and wouldn't be homogeneous.) Writing $\alpha=a+i d$ then yields

$$f(z)=(a+i d)(x+i y)^3=(ax^3-3dx^2y-3axy^2+dy^3)+i(dx^3+3ax^2y-3dxy^2-ay^3)$$ in agreement with $u(x,y)$, $v(x,y)$ as identified in the OP.