An efficient pseudospectral method for performing fully-coupled six-dimensional bound state dynamics calculations is presented. A Lanczos-based iterative diagonalization scheme produces the energy levels in increasing energies. This scheme, which requires repetitively acting the Hamiltonian operator on a vector, circumvents the problem of constructing the full matrix. This permits the use of ultralarge molecular basis sets in order to fully converge the calculations. The Lanczos scheme was conducted in a symmetry adapted six-dimensional spectral representation. The Hamiltonian operator has been split into only four different terms, each being Hermitian and symmetry-adapted. The potential term is evaluated by a pseudospectral scheme of Gaussian accuracy, which guarantees the variational principle. Spectroscopic levels are computed with this method for one ammonia potential, and compared to experimental results. The results presented below are a direct application of our vector formulation. The latter has shown to be particularly well adapted to the split pseudospectral approach for it yields a compact and symmetry-adapted Hamiltonian.