Fundamentally, here are the building blocks of what you are asking for:

* Consider solving the problem
$$
  -\Delta u = f
$$
in a domain $\Omega$ with boundary values $u=g$ on $\partial\Omega$.

* **Theorem 1:** For any $\Omega,f,g$, the solution can be written as
$$
  u(\mathbf x) = \sum_{i=1}^\infty A_i \varphi_i(\mathbf x)
$$
where the functions $\varphi_i$ are the *eigenfunctions* of the Laplace operator and satisfy
$$
  -\Delta \varphi_i = \lambda_i \varphi_i,
$$
and the $A_i$ are expansion coefficients that can be computed via $A_i = \frac{1}{\lambda_i}\int_\Omega \varphi_i(\mathbf x) f(\mathbf x)\, dx$ plus some terms that come from the boundary values $g$. (The theorem is true because the eigenfunctions of the Laplace operator form a *complete basis* of $L_2(\Omega)$ and consequently of the solution space in which $u$ lies.)

* **Theorem 2:** For a given $\Omega,f,g$, the coefficients $A_i$ must decay, that is: $A_i\rightarrow 0$. As a consequence, you can *approximate* the solution by truncating the sum as
$$
  u(\mathbf x) \approx \sum_{i=1}^N A_i \varphi_i(\mathbf x).
$$
(This theorem is true because the norm of $u$ must be finite, and that can only be the case if the sum $\sum_{i=1}^\infty A_i$ converges -- which it only does if the coefficients decay sufficiently fast.)

* **Theorem 3:** The eigenfunctions $\varphi_i$ can be approximated numerically, for example using the finite element method, by functions $\varphi_i^h$. Then the solution of the original problem can be approximated as
$$
  u(\mathbf x) \approx u^h(\mathbf x) = \sum_{i=1}^N A_i \varphi_i^h(\mathbf x).
$$

I will note, however, that using this expansion into eigenfunctions *is not an efficient way to find an approximate solution of the Laplace equation*.