I was playing around with Laplace's equation:
$$\frac{\partial^2 u}{\partial x ^2}+\frac{\partial^2 u}{\partial y ^2}=0$$ It occured to me that it can be written as:
$$\bigg(\frac{\partial}{\partial x}-\frac{1}{i}\frac{\partial}{\partial y}\bigg)\bigg(\frac{\partial}{\partial x}+\frac{1}{i}\frac{\partial}{\partial y}\bigg)u=0$$
That means that the most general solution to the equation is:
$$u=f(x+iy)+g(x-iy)$$
That is:
$$u=f(z)+g(\bar z)$$
Suppose we want to solve the equation in a finite domain of $\mathbb{R}^2$ with Dirichlet boundary conditions. Isn't it possible to reduce the problem to an integral using Cauchy's formula of complex analysis?
Yes, the Laplace equation in two dimensions is often written as $u_{z\bar z}=0$, which simplifies some computation. And yes, a harmonic function can be locally represented as the sum of a holomorphic and antiholomorphic functions. In a simply-connected domain this representation is also global; but in general multi-valued functions may enter the picture. For example, consider $u(z)=\log |z|$, a harmonic function in the annulus $\Omega = \{1<|z|<2\}$. We can write $$u(z)=\frac12\log z+\frac12\log\bar z$$ but both terms are multi-valued in $\Omega$. The representation is still useful.
The second part of your question is about boundary value problems. Sure, holomorphic functions can be recovered from boundary values by means of the Cauchy formula, using the kernel $\frac{1}{z-\zeta}$. This looks easier than recovering harmonic functions with the Poisson kernel, which depends on the domain and can rarely be found explicitly.
However, there is a catch. Cauchy's formula requires the knowledge of both real and imaginary parts of function, which is an overdetermined boundary condition for holomorphic functions. In terms of harmonic functions, this means having the boundary values for both $u$ and its conjugate $\tilde u$. But we don't have the latter with the Dirichlet problem for harmonic functions, so the recovery is more difficult. Formal manipulations with $f(z)+g(\bar z)$ don't give us more information: we don't know the boundary values of either $f$ or $g$.
Put another way: the Cauchy kernel is simple because it recovers a rigid object (a holomorphic function), while the Poisson kernel has to take into account the geometry of the domain which affects the behaviour of a harmonic function.