Suppose I want to solve the Schrödinger equation numerically for some potential V(x). The easiest way to do so, is to discretize it on a grid of finite length, and apply a finite difference scheme to approximate the second order derivative. Doing so yields an eigenvalue equation on matrix form for the wavefunctions and their corresponding energies, which may then be found by diagonalization. The above method was the way I always solved the Schrödinger equation, when numerical work was needed. However, now I am faced with a problem, where I need to impose a boundary condition on the boundary of my interval. In the method above you implicitly assume that the wavefunctions outside the domain you are looking at, but for my current problem this will no longer work. Is there a way to adapt the finite difference method above to handle the case with a non-zero boundary condition? Worse even, my boundary condition differs for the different eigenmodes. I.e. the groundstate has one value at the boundary, the first excited another and so on. Could this also be incorporated easily?