# Numerical solution to the infinite well problem

I've used the following code to implement it

import numpy as np
from scipy.linalg import eig
import matplotlib.pyplot as plt

W = 10 #width of well
n = 1000 #number of points excluding 0,W
hbar = 1;m=1 #atomic units
x = np.linspace(1,W-1,num=n) #exclusion of 0,W
h = x[1]-x[0]#step size

#Construction of the hamiltonian matrix
H = np.zeros([n,n],dtype=complex)
for i in range(n):
H[i][i] = -2
if(i-1>=0):
H[i][i-1] = 1
if(i+1<n):
H[i][i+1] = 1
H = -1*(hbar**2)/(2*m*(h**2))*H
E,V = eig(H)

plt.plot(x,np.abs(V[:][0]),'.')
plt.plot(x,np.abs(V[:][1]),'.')
plt.show()



Result:

I'm wondering why exactly this code is not optimal enough to reproduce results close to the analytical solutions, and what the source of error is.
Edit:
Output for N=10,000

However,it took about 45 min for it to run

• What happens when you increase the discretization (number of points)? If it becomes any better, consider using Richardson interpolation here. Jul 5 at 15:06
• I think that is a good idea that you write the differential equation and finite differences used. That will make easier to understand the problem at hand for more people. Jul 5 at 15:50

I think that the main problem might be with the solver you are using. The Hamiltonian (matrix) in this case is Hermitian, it is even symmetric since it is purely real. You could use eigh instead of eig to take advantage of this.

Furthermore, you are not removing only the first and last points but intervals of size 1 at each end. Following, I show you a snippet with this changes.

import numpy as np
from scipy.linalg import eigh
import matplotlib.pyplot as plt

W = 10 # width of well
n = 1000 # number of points
hbar = 1
m = 1
x = np.linspace(0, W, num=n+2)
h = x[1] - x[0] # step size

# Construction of the hamiltonian matrix
H = np.zeros([n, n])
for i in range(n):
H[i, i] = -2
if(i-1>=0):
H[i, i-1] = 1
if(i+1<n):
H[i, i+1] = 1
H = -1*(hbar**2)/(2*m*(h**2))*H

# Solution
E, V = eigh(H, eigvals=(0, 10))

# Visualization
eigfun = np.zeros((n+2))
eigfun[1:-1] = V[:, 0]
plt.plot(x, np.abs(eigfun)**2)
eigfun[1:-1] = V[:, 1]
plt.plot(x, np.abs(eigfun)**2)
plt.show()


Also, you are visualizing $$|\psi|$$, but I think that it is better to visualize $$|\psi|^2$$.

Finally, you could take advantage of the sparsity of the matrix. See this answer.