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
        H[i][i-1] = 1
        H[i][i+1] = 1
H = -1*(hbar**2)/(2*m*(h**2))*H 
E,V = eig(H)


Result: enter image description here

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.
Output for N=10,000

However,it took about 45 min for it to run

  • $\begingroup$ What happens when you increase the discretization (number of points)? If it becomes any better, consider using Richardson interpolation here. $\endgroup$
    – Paddy
    Jul 5, 2021 at 15:06
  • $\begingroup$ 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. $\endgroup$
    – nicoguaro
    Jul 5, 2021 at 15:50

1 Answer 1


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
        H[i, i-1] = 1
        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)

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

enter image description here

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.