In the present paper, the authors consider the Schrödinger operator H with the Coulomb potential defined in R 3m , where m is a positive integer. Both bounded domain approximations to multielectron systems and finite element approximations to the helium system are analyzed. The spectrum of H becomes completely discrete when confined to bounded domains. The error estimate of the bounded domain approximation to the discrete spectrum of H is obtained. Since numerical solution is difficult for a higher-dimensional problem of dimension more than three, the finite element analyses in this paper are restricted to the S-state of the helium atom. The authors transform the six-dimensional Schrödinger equation of the helium S-state into a three-dimensional form. Optimal error estimates for the finite element approximation to the three-dimensional equation, for all eigenvalues and eigenfunctions of the three-dimensional equation, are obtained by means of local regularization. Numerical results are shown in the last section. 1. Introduction. The multielectron Coulomb problem in quantum mechanics cannot be solved in a finite form. Nevertheless it challenges and stimulates many mathematicians and physicists to devote themselves to developing efficient methods for solving the system. Several successful approximation techniques in quantum physics/chemistry have been developed for this problem. They include the Hartree-Fock method [15], the finite difference method [19], [35], the correlation-function hyperspherical-harmonic method [18], [24], and various variational approximations. For the Hartree-Fock method, every electron is considered independently to be in a central electric field formed by the nucleus and other electrons. The finite difference method needs a rectangular domain in R N and uniform grids. The double and triple basis set methods (which are variational methods indeed) are very powerful for the eigenvalue problem of the helium atom. Kono and Hattori (see [21], [22]) used two sets of basis functions r i 1 r j 2 r k 12 e −ξr1−ηr2 A ("ξ terms") and r i 1 r j 2 r k 12 e −ζ(r1+r2) A ("ζ terms") to calculate the energy levels for the S, P , and D states of the helium atom. (A is an appropriate angular factor.) The former set of functions is expected to describe the whole wave function roughly, while the latter is expected to describe the short-and middle-range correlation effects. Their calculations yield 9-10 significant digits for S states. Kleindienst, Lüchow, and Merckens [20] and Drake and Yan [12] applied the double basis set method to S-states of helium. Their basis functions are r i 1 r j 2 r k 12 e −ξ1r1−η1r2 and r i 1 r j 2 r k 12 e −ξ2r1−η2r2 . Drake and Yan employed truncations to ensure numerical stability and convergence. By complete optimization of the exponential scale factors ξ 1 ,