We introduce an accurate and efficient algebraic technique for the computation of the vibrational spectra of triatomic molecules, of both linear and bent equilibrium geometry. The full three-dimensional potential energy surface (PES), which can be based on entirely ab initio data, is parametrized as a product Morse-cosine expansion, expressed in bond angle internal coordinates, and includes explicit interactions among the local modes. We describe the stretching degrees of freedom in the framework of a Morse-type expansion on a suitable algebraic basis, which provides exact analytical expressions for the elements of a sparse Hamiltonian matrix. Likewise, we use a cosine power expansion on a spherical harmonics basis for the bending degree of freedom. The resulting matrix representation in the product space is very sparse, and vibrational levels and eigenfunctions can be obtained by efficient diagonalization techniques. We apply this method to carbonyl sulfide, hydrogen cyanide, water, and nitrogen dioxide. When we base our calculations on high-quality PESs tuned to the experimental data, the computed spectra are in very good agreement with the observed band origins.