I am right now working on a script that solves the Schrodinger equation numerically for arbitrary potentials using the finite difference method. The idea is that I diagonalize the Hamiltonian with elements: H(i,i+1)=1/dx^2 * constants H(i,i-1)=1/dx^2 * constants H(i,i) = -2/dx^2 * constants and zero elsewhere. Now what I am not so sure about is what boundary conditions this procedure assumes? I think it assumes automatically that the wave functions tend to zero at the boundaries of your grid. In any case the script works for 1d, but I am now trying to make it work for 2d so I can solve the circular well problem. However, it has not worked so far: I have used that the radial equation (polar coordinates) is given as that shown here http://www.physics.csbsju.edu/QM/square.08.html which has led me to a hamiltonian of the form (dx is mesh size) H(i,i) = -2/dx^2 - 1/r(i)^2 H(i,i+1)=1/dx^2 + 1/(dx*r(i)) H(i,i-1)=1/dx^2 - 1/(dx*r(i)) But the solutions I am getting oscillate more and more and don't look anything like the analytic solution. Can anyone see what the problem is? I suspect it may be the fact that my method forces the wave function to go to zero at the boundaries, which does not make sense for r=0. If this is indeed the problem, what could I do to fix it?