We present a nonlinear eigenvalue solver enabling the calculation of bound solutions of the Schrodinger equation in a system with contacts. We discuss how the imposition of contacts leads to a nonlinear eigenvalue problem and discuss the numerics for a one- and two-dimensional potential. We reformulate the problem so that the eigenvalue problem can be efficiently solved by the recently proposal rational Krylov method for nonlinear eigenvalue problems, known as NLEIGS. In order to improve the convergence of the method, we propose a holomorphic extension such that we can easily deal with the branch points introduced by a square root. We use our method to determine the bound states of the one-dimensional Poschl-Teller potential, a two-dimensional potential describing a particle in a canyon and the multi-band Hamiltonian of a topological insulator.