We present a stable scheme to calculate continuum and bound electronic states in the vicinity of a surface of a semi-infinite crystal. The method is based on the Numerov integration algorithm and uses singular value decomposition to control the exponentially growing contributions. The method has been successfully tested on the Li (110) surface with and without adsorbed H. For this model system, we are able to locate the energies of H-induced surface states relative to the corresponding energies of bulk continuuum states. Results encourage further development.