  With the aim of describing bound and continuum states for diatomic molecules,
we develop and implement a spectral method that makes use of Generalized
Sturmian Functions (GSF) in prolate spheroidal coordinates. In order to master
all computational issues, we apply here the method to one--electron molecular
ions and compare it with benchmark data for both ground and excited states. We
actually propose two different computational schemes to solve the two coupled
differential equations.
  The first one is an iterative 1d procedure in which one solves alternately
the angular and the radial equations, the latter yielding the state energy. The
second, named direct $2d$ method, consists in representing the Hamiltonian
matrix in a two--dimensional GSF basis set, and its further diagonalization.
Both spectral schemes are timewise computationally efficient since the basis
elements are such that no derivatives have to be calculated numerically.
Moreover, very accurate results are obtained with minimal basis sets. This is
related on one side to the use of the natural coordinate system and, on the
other, to the intrinsic good property of all GSF basis elements that are
constructed as to obey appropriate physical boundary conditions. The present
implementation for bound states paves the way for the study of continuum states
involved in ionization of one or two--electron diatomic targets.
