Suppose I am given some 1D Hamiltonian: H = ħ2/2m d2/dx2 + V(x) (1) Which I want to solve on the interval [0,L]. I think most of you are familiar with the standard approach of discretizing the interval [0,L] in N pieces and using the finite difference formulas for V and the second derivative in (1), which can then be formulated as a matrix equation, which may be diagonalized for the eigenvector and eigenvalues. Now for all this to work one has to assume that the wave-function goes to zero outside the interval [0,L], which follows if one makes the effort of writing up the finite difference expressions. My question is: Is there a way to enforce another boundary condition with this method? I am solving a problem, where it would be beneficial to enforce the wave function to take a non-zero value on the boundaries of the interval [0,L]. Is this possible with the standard finite difference method or should I look at more advanced methods?