In this paper, we address the problem of determining and efficiently computing an approximation to the eigenvalues of the negative Laplacian -Delta on a general domain Omega subset of R-2 subject to homogeneous Dirichlet or Neumann boundary conditions. The basic idea is to look for eigenfunctions as the superposition of generalized eigenfunctions of the corresponding free space operator, in the spirit of the classical method of particular solutions (MPS). The main novelties of the proposed approach are the possibility of targeting each eigenvalue independently without the need for extensive scanning of the positive real axis and the use of small matrices. This is made possible by iterative inclusion of more basis functions in the expansions and a projection idea that transforms the minimization problem associated with MPS and its variants into a relatively simple zero-finding problem, even for expansions with very Jew basis functions.