In this paper, we study natural boundary reduction for Laplace equation with Dirichlet or Neumann boundary condition in a three-dimensional unbounded domain, which is the outside domain of a prolate spheroid. We express the Poisson integral formula and natural integral operator in a series form explicitly. Thus the original problem is reduced to a boundary integral equation on a prolate spheroid. The variational formula for the reduced problem and its well-posedness are discussed. Boundary element approximation for the variational problem and its error estimates, which have relation to the mesh size and the terms after the series is truncated, are also presented. Two numerical examples are presented to demonstrate the effectiveness and error estimates of this method.