We present a contracted basis-iterative method for calculating numerically exact vibrational energy levels of methane (a 9D calculation). The basis functions we use are products of eigenfunctions of bend and stretch Hamiltonians obtained by freezing coordinates at equilibrium. The basis functions represent the desired wavefunctions well, yet are simple enough that matrix-vector products may be evaluated efficiently. We use Radau polyspherical coordinates. The bend functions are computed in a nondirect product finite basis representation [J. Chem. Phys. 118, 6956 (2003)] and the stretch functions are computed in a product potential optimized discrete variable (PODVR) basis. The memory required to store the bend basis is reduced by a factor of ten by storing it on a compacted grid. The stretch basis is optimized by discarding PODVR functions with high potential energies. The size of the primitive basis is 33 billion. The size of the product contracted basis is six orders of magnitude smaller. Parity symmetry and exchange symmetry between two of the H atoms are employed in the final product contracted basis. A large number of vibrational levels are well converged. These include almost all states up to 8000 cm(-1) and some higher local mode stretch bands. (C) 2003 American Institute of Physics.