  Let $\Omega$ be an open, simply connected, and bounded region in
$\mathbb{R}^{d}$, $d\geq2$, and assume its boundary $\partial\Omega$ is smooth.
Consider solving an elliptic partial differential equation $-\Delta u+\gamma
u=f$ over $\Omega$ with a Neumann boundary condition. The problem is converted
to an equivalent elliptic problem over the unit ball $B$, and then a spectral
Galerkin method is used to create a convergent sequence of multivariate
polynomials $u_{n}$ of degree $\leq n$ that is convergent to $u$. The
transformation from $\Omega$ to $B$ requires a special analytical calculation
for its implementation. With sufficiently smooth problem parameters, the method
is shown to be rapidly convergent. For $u\in C^{\infty}(\overline{\Omega}) $
and assuming $\partial\Omega$ is a $C^{\infty}$ boundary, the convergence of
$\Vert u-u_{n}\Vert_{H^{1}}$ to zero is faster than any power of $1/n$.
Numerical examples in $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$ show experimentally
an exponential rate of convergence.
