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.