Factorization methods such as the Hamiltonian hierarchy have been useful to find eigenfunctions for Schrödinger equations, in particular, for potentials that are partially or approximately solvable. In this paper, an alternative approach is proposed to study excited states via the variational method. The trial functions are built from the exact or approximate superpotential for the ground state combined with the Gram–Schmidt process to ensure orthogonalization between the functions. The results found variationally for one dimensional potentials are compared with previous results from the literature. The energy eigenvalues obtained agree with previous ones and, for most of the results, the percentage difference between the proposed approach and others in the literature is less than 0.1%. The method introduced is an effective and intuitive approach to determine trial wave functions for the excited states. This approach can be useful in studying the Schrödinger equation and related problems which can be mapped onto a Schrödinger type-equation as, for example, the Fokker–Planck equation.