Woods-Saxon energy levels

I'm trying to calculate the bound state energy levels for a Woods-Saxon potential, given by:

$$V(r) = -\frac{V_{0}}{1+\exp(\frac{r-R}{a})}$$

I'm using Numerov's matrix method, which led me to:

$$H\psi = E\psi, \hspace{2cm} H= -\frac{\hbar^{2}}{2m}B^{-1}A + V$$

with:

$$A = \frac{1}{d^{2}}\begin{pmatrix} -2 & 1 & 0 & 0 & 0 & ...\\ 1 & -2 & 1 & 0 & 0 & ...\\ 0 & 1 & -2 & 1 & 0 & ... \\ 0 & 0 & 1 & -2 & 1 & ... \\ ... & ... & ... & ... & ... & ... \\ \end{pmatrix}$$

$$B = \frac{1}{12}\begin{pmatrix} 10 & 1 & 0 & 0 & 0 & ...\\ 1 & 10 & 1 & 0 & 0 & ...\\ 0 & 1 & 10 & 1 & 0 & ... \\ 0 & 0 & 1 & 10 & 1 & ... \\ ... & ... & ... & ... & ... & ... \\ \end{pmatrix}$$

$$V = \begin{pmatrix} V_{1} & 0 & 0 & 0 & 0 & ...\\ 0 & V_{2} & 0 & 0 & 0 & ...\\ 0 & 0 & V_{3} & 0 & 0 & ... \\ 0 & 0 & 0 & V_{4} & 0 & ... \\ ... & ... & ... & ... & ... & ... \\ \end{pmatrix}$$

$$\psi = \begin{pmatrix} \psi_{1}\\ \psi_{2}\\ \psi_{3}\\ \psi_{4}\\ ...\\ \end{pmatrix}$$

as made in Bhat et al (2015).

However, this is an eigenvalue equation, and I get $$n$$ eigenvalues, one for each point in space $$(r_{n} = n\cdot dr)$$, and I don't know what to do with them. How can I know which ones are the allowed energy states? I think I should set some boundary conditions for the wave function, but I'm not sure how to do it. Here's most of my Python code, the rest is just plots and stuff.

import numpy as np
import scipy as sp
import scipy.linalg
from scipy.sparse import diags

def WS(r, V0, R, a):
return -V0/(1.0+sp.exp((r-R)/a))
hc = 197
m = 939.57   # MeV/c2
V0 = 50   # MeV
a = 0.6   # fm
R = 5.0   # fm

dr = 0.1
r0 = 0.0
n = 100
r = np.zeros(n)
for i in range(n):
r[i] = r0 + i*dr

A = diags([np.full(n,-2), np.full(n-1,1), np.full(n-1,1)], [0, -1, 1]).toarray()/(dr**2)
B = (1.0/12.0)*diags([np.full(n,10), np.full(n-1,1), np.full(n-1,1)], [0, -1, 1]).toarray()
V = np.zeros((n,n))
for i in range(n):
V[i,i] = WS(r[i], V0, R, a)
H = (-hc**2/(2*m))*np.dot(np.linalg.inv(B), A) + V

E = np.linalg.eig(H)[0]
psi = scipy.linalg.solve(H, E)


So the big question is: how can I get the bound-state allowed energies? What am I missing, or what did I wrong?

• I'm not reading all that code but are you treating the energy level as a variable? This looks like a numerical method for solving an ODE or PDE and you either need the constant E or you need it to fall out. The standard procedure is to make E a variable and add dE/dx = 0, dE/dt = 0, etc. This works well with the relaxation method where you can start with a guess involving exact solutions to another problem. If this is not the issue someone else may be able to help. – ggcg Dec 27 '18 at 17:32
• Yeah, the eigenvalues of the discrete matrix are not necessarily related to the eigenvalues of the differential operator. That is your first mistake. – ggcg Dec 27 '18 at 17:33