Excitonic effects in optical spectra and electron-hole pair excitations are described by solutions of the Bethe-Salpeter equation (BSE) that accounts for the Coulomb interaction of excited electron-hole pairs. Although for the computation of excitoni
c optical spectra in an extended frequency range efficient methods are available, the determination and analysis of individual exciton states still requires the diagonalization of the electron-hole Hamiltonian $hat{H}$. We present a numerically efficient approach for the calculation of exciton states with quadratically scaling complexity, which significantly diminishes the computational costs compared to the commonly used cubically scaling direct-diagonalization schemes. The accuracy and performance of this approach is demonstrated by solving the BSE numerically for the Wannier-Mott two-band model in {bf k} space and the semiconductors MgO and InN. For the convergence with respect to the $vk$-point sampling a general trend is identified, which can be used to extrapolate converged results for the binding energies of the lowest bound states.
