ground state from the Schroedinger equation with a central potential what happens to the origin

Question

I have code that attempts to implement a solution to the Schrödinger equation where there is a central potential (more or less im thinking of hydrogen), in 1-D using the numerov method to construct the wave function. It does not take care of the states (N=1,..) as my goal is to understand the link between the theory and the software. Here are some results for some energies, whether or not they are eigenvalues, im not sure:

[6.054999999999915, 6.214999999999912, 6.534999999999905]

and a graph. The algorithm im trying to implement is noted here in section 6 PDF, but my potential is different. My problem is, the central potential has a 1/r at the origin, how would I deal with that? (I did a sidestep in the code by adding an arbitrary 0.00001 to eliminate 0).

There is a more general code in another question seen here, but I do not follow the shift from equation to 3d matrix.

Here is my code:

import numpy as np
import scipy
import matplotlib.pyplot as plt 

h=np.float64(1)
m=60
c=m/2
estep = 1.0/(100)
delta = 0.0002
x=np.linspace(-1000,1000,num=m)
x=(x - x[m/2]+0.00001 )
v =  1.0/(x**2)
sig=None
E = []

def get_left(wf,k):
    for i in range(2,int(c)):
        wf[i] = 2*(( 1- 1/12* h**2* k[i-2]**2)*wf[i-2] - (1 + 1/12 *h**2* k[i-1]**2)*wf[i-1])/(1 + 1/12* h**2* k[i]**2)
    wf[0:c] = wf[0:c]/(np.sqrt(np.sum(wf[0:c]**2)))
    return wf

def get_right(wf,k):
    mi = m-1
    for i in range(2,int(c)):
        wf[mi-i] = 2*(( 1- 1/12* h**2* k[mi-i+2]**2)*wf[mi-i+2] - (1 + 1/12 *h**2* k[mi-i+1]**2)*wf[mi-i+1])/(1 + 1/12* h**2*k[mi-i]**2)
    wf[c:mi] = wf[c:mi]/(np.sqrt(np.sum(wf[c:mi]**2)))
    return wf

def root_search(epsilonu, epsilonl):
    en = None
    it =0
    while it < 10000:
        phil =  np.zeros(m)
        phir =  np.zeros(m)
        phil[0] =0
        phil[1] = 0.001
        phir[m-1]=0
        phir[m-2]=0.001
        k = (-v + (epsilonu+epsilonl)/2)
        phil = get_left(phil,k)
        phir = get_right(phir,k)
        diffl = np.diff(phil)
        diffr = np.diff(phir)
        erri = diffl[m/2]- diffr[m/2]
        if erri < delta:
            plt.plot(x, (phil+phir)**2 )
            plt.savefig( str('sho'+'.png' ))
            return (epsilonu+epsilonl)/2
        if erri < 0 : 
            epsilonl = (epsilonu+epsilonl)/2
        if erri > 0:
            epsilonu =epsilonl 
            epsilonl = (epsilonu+epsilonl)/2
        it+=1
    return None

def start():
    epsilonl=0.01
    eps = epsilonl
    epsilonu=None
    iteration=0
    sig = None
    while iteration < 100000:
        phil =  np.zeros(m)
        phir =  np.zeros(m)
        phil[0] =0
        phil[1] = 0.001
        phir[m-1]=0
        phir[m-2]=0.001
        k = (-v + eps)
        phil = get_left(phil,k)
        phir = get_right(phir,k)
        diffl = np.diff(phil)
        diffr = np.diff(phir)
        err = diffl[m/2]- diffr[m/2]
        ssig = 1 if err>0 else -1
        if sig:
            if sig != ssig:
                epsilonu= eps 
                E.append(root_search(epsilonu, epsilonl))
                print(E)
        sig = ssig
        epsilonl = eps
        eps = eps+estep
        iteration+=1

if __name__ == "__main__":
    start()

The theory: Im attempting to solve the 1-D Schroedinger equation, $$ \left( -\frac{h^2}{2m}\frac{d^2}{dx^2} + V\left(x\right)\right) \phi\left(x\right) = \epsilon \phi\left( x\right) $$ with V as a central potential (coulombic force, inverse square relation), the final equation looks like (with proper choice of units to give q^2=1, h-1..., energy in hartree) $$ \frac{d^2}{dx^2}\phi\left(x\right) + 2 \left(\epsilon -\frac{1}{x^2} \right) \phi\left(x\right) =0 $$

Answer 1

As mentioned by @Kirill, the differential equation for the hydrogen is slightly different. After the separation of variables, you end up with

$$\left[ - \frac{\hbar^2}{2\mu} \left({1 \over r^2}{\partial \over \partial r}\left(r^2 {\partial R(r)\over \partial r}\right) - {l(l+1)R(r)\over r^2} \right) + V(r)R(r) \right]= E R(r),$$

You can rewrite the equation as

$$\left[ -\frac{1}{2} \frac{d^2}{dy^2} + \frac{1}{2} \frac{l(l+1)}{y^2} - \frac{1}{y}\right] u_l = W u_l$$

The next code (based on my previous post) computes the eigenvalues/eigenvectors for this equation. See that I added a small value to the denominators, so they don't vanish at zero.

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

L = 50*np.pi
N = 1001
xa = 0
xb = L
x = np.linspace(xa, xb, N)
dx = x[1] - x[0]

l = 1
k = 0
T = -0.5*diags([-2., 1., 1.], [0, -1, 1], shape=(N, N))/dx**2
U_vec = 0.5*l*(l + 1)/(x**2 + 1e-6) - 1/(np.abs(x) + 1e-6)
U = diags([U_vec], [0])

H = T + U

vals, vecs = eigsh(H, which='SA')

print(np.round(vals, 6))
print(np.round([-1/(2*n**2) for n in range(k + l + 1, k + l + 7)], 6))

for k in range(5):
    vec = vecs[:, k]
    mag = np.sqrt(np.dot(vecs[:, k],vecs[:, k]))
    vec = vec/mag
    plt.plot(x, vec, label=r"$n=%i$"% (k+1))

plt.xlabel(r"$x$")
plt.ylabel(r"$\psi(x)$")
plt.xlim(xa, xb)
plt.legend()
plt.savefig("eigenvectors.png", dpi=600) 
plt.show()

That gives as solution:

[-0.125106 -0.0556   -0.031272 -0.020012 -0.013896 -0.010209]
[-0.125    -0.055556 -0.03125  -0.02     -0.013889 -0.010204]

For your repulsive potential (different sign than the Laplacian), I used the same code

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

L = 50*np.pi
N = 1001
xa = 0
xb = L
x = np.linspace(xa, xb, N)
dx = x[1] - x[0]

T = -0.5*diags([-2., 1., 1.], [0, -1, 1], shape=(N, N))/dx**2
U_vec = 1/(x**2 + 1e-6)
U = diags([U_vec], [0])

H = T + U

vals, vecs = eigsh(H, which='SA')

print(np.round(vals, 6))

for k in range(5):
    vec = vecs[:, k]
    mag = np.sqrt(np.dot(vecs[:, k],vecs[:, k]))
    vec = vec/mag
    plt.plot(x, vec, label=r"$n=%i$"% (k+1))

plt.xlabel(r"$x$")
plt.ylabel(r"$\psi(x)$")
plt.xlim(xa, xb)
plt.legend()
plt.savefig("eigenvecs.png", dpi=600)  
plt.show()

with answers

[ 0.000408  0.001207  0.002405  0.004001  0.005997  0.008392]

Although, I am not sure what that means. Because it is not letting the electron to be in the origin pushing it away, and the domain should be really unbounded. It is similar to having a square potential with a round edge at the left. For the more realistic case, you can use nonuniform finite differences and increase them exponentially or absorbing boundary conditions.