We present a spectral method for solving elliptic equations which arise in general relativity, namely three-dimensional scalar Poisson equations, as well as generalized vectorial Poisson equations of the type $��\vec{N} + ��\vec{\nabla}(\vec{\nabla}\cdot \vec{N}) = \vec{S}$ with $��\not= -1$. The source can extend in all the Euclidean space ${\bf R}^3$, provided it decays at least as $r^{-3}$. A multi-domain approach is used, along with spherical coordinates $(r,��,��)$. In each domain, Chebyshev polynomials (in $r$ or $1/r$) and spherical harmonics (in $��$ and $��$) expansions are used. If the source decays as $r^{-k}$ the error of the numerical solution is shown to decrease at least as $N^{-2(k-2)}$, where $N$ is the number of Chebyshev coefficients. The error is even evanescent, i.e. decreases as $\exp(-N)$, if the source does not contain any spherical harmonics of index $l\geq k -3$ (scalar case) or $l\geq k-5$ (vectorial case).

Minor revisions. 31 pages, 13 figures, accepted for publication J. Comp. Phys