$u(x,y)=e^{x^2-y^2}(e^y \cos(x-2xy)+ e^{-y} \cos(x+2xy))$
Solving the Cauchy-Riemann equations for this is not practical. Alternatively, I could try to express $u(x,y)$ in the form of the real or complex part of a (analytic) function in $z$. In doing that, I either observe it right away, or put in $x=\frac{z+\bar z}{2}$ and $y=\frac{z-\bar z}{2i}$. But that turned really messy and the $z$'s and $\bar z$'s in the coefficients of the exponential functions do not reduce too give the desired form. I would like to ask what else can I do here to find the harmonic conjugate?
Given that $u$ is harmonic, you can simply replace $x$ with $z$ and $y$ with $0$ to get an analytic function which agrees with $u$ on the real axis, and hence everywhere: $$ f(z) = 2 e^{z^2} \cos z . $$ Now just take the imaginary part of $f(x+iy)$.