A variational formulation finite element method is developed for calculation of vibrational wave functions in a domain spanned by close-coupled, or Jacobi, coordinates R and gamma. C-1 tensor-product basis functions, which allow straightforward separation of kinetic and overlap integrals into products of one-dimensional integrals, are used. Furthermore, representation of the potential energy surface in terms of the same tensor-product basis functions used to represent the wave functions allows the potential energy integrals to also be written as a sum of products of one-dimensional integrals. Factorization of the integrals, together with expression of one-dimensional integrals in analytic or rapidly convergent power series form, reduces the computational effort of calculation of all matrix elements to a small, and arguably insignificant, level. It is shown that the theoretical error in eigenvalue, i.e., O(h(6)) for bicubic Hermite functions, is achieved for a number of rare gas van der Waals triatomics for which surfaces have been previously published. We also present illustrative calculations on NeHCl and (2)A(') and (2)A(') NeHCl+, which have not been previously studied, for surfaces calculated at the CCSD(T)/cc-pVTZ level.