# Solve the stationary Schrödinger equation with the finite difference method [closed]

I am trying to solve the stationary Schrödinger equation for a double-level well potential through the finite difference method.

Here is the shape of the potential I would like to solve it for where the second walls are infinite.

The Schrödinger equation for this problem is the following

$$-\frac{\hbar^2}{2m} \partial_{xx}\psi(x) = [E-V(x)]\psi(x)$$

which can be turned into

$$\partial_{xx}\psi(x) = \frac{2m}{\hbar^2}[V(x)-E]\psi(x)$$

we can discretise the wave function and approximate the second derivative of $$\psi_i$$ as

$$\partial_{xx}\psi_i \sim \frac{\psi_{i-1} - 2\psi_i + \psi_{i+1}}{h^2}$$

Then the Schrödinger equation becomes

$$\frac{\psi_{i-1} - 2\psi_i + \psi_{i+1}}{h^2} = \frac{2m}{\hbar^2}[V(x_i)-E]\psi(x)$$

which can be put as an eigenvalue problem as this

$$\begin{pmatrix} -2 & 1 & 0 & ... & 0\\ 1 & -2 & 1 & ... & 0\\ 0 & .&.&. & 0\\ 0 & ... & 1 & -2 & 1 \end{pmatrix} \cdot \begin{pmatrix} \psi_1\\ \psi_2\\ .\\ .\\ .\\ \psi_N \end{pmatrix} = \frac{2m}{\hbar^2} \begin{pmatrix} V_1 - E & 0 & ... & 0\\ 0 & V_2 - E & ... & 0\\ 0 & 0 & ... & 0\\ 0 & 0 & ... & V_N - E \end{pmatrix} \cdot \begin{pmatrix} \psi_1\\ \psi_2\\ .\\ .\\ .\\ \psi_N \end{pmatrix}$$

Here is my question:

If I didn't have $$V_i$$ I would be using

scipy.linalg.eigh()

to solve this problem for $$E$$, however $$E$$ is now inside the second matrix. Is there a way to solve this equation finding the values of $$E$$ and relative sets of eigenvectors with Python?

I implemented it as suggested in the comments as

$$\begin{pmatrix} -2-\frac{2m}{\hbar^2}V_1 & 1 & 0 & ... & 0\\ 1 & -2-\frac{2m}{\hbar^2}V_2 & 1 & ... & 0\\ 0 & .&.&. & 0\\ 0 & 0 &... & 1 & -2-\frac{2m}{\hbar^2}V_N \end{pmatrix} \cdot \begin{pmatrix} \psi_1\\ \psi_2\\ .\\ .\\ .\\ \psi_N \end{pmatrix} = -\frac{2m}{\hbar^2}E \begin{pmatrix} \psi_1\\ \psi_2\\ .\\ .\\ .\\ \psi_N \end{pmatrix}$$

and used scipy.linalg.eigh().

Here is my code

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

# number of discrete bins in which the domain is divided
N = 200

# create discretised arrays for the x axis and the value of the wave funciton
xaxis = np.linspace(-1, 1, N)
Vpot = np.linspace(-1, 1, N)

# factor with energy hbar and mass
m = 1
hbar = 1

# define the potential
def V(x):

if np.abs(x)<0.5:
val = 0.
else:
val = 4.

return val

# fill an array for the potential
for i in range(len(xaxis)):
Vpot[i] = V(xaxis[i])

# initialize NxN matrix
Hmat = [[0 for x in range(N)] for y in range(N)]

# fill Hmat like the following
#   -2-2m/h^2*V1 1 0 0 ...0
#   1 -2-2m/h^2*V2 1 0 ...0
#   0 1 -2-2m/h^2*V3 1 ...0
#   .....................0
#   .........1 -2-2m/h^2*VN
for row in range(N):
for elem in range(N):
if row == elem: Hmat[row][elem] = -2 -(2*m)/(hbar**2)*Vpot[elem]
if np.abs(row-elem) == 1: Hmat[row][elem] = 1

# now I get the eigenvalue with w[i] and corresponding eigenvector v[i]
energies, psi = eigh(Hmat, b=None, eigvals_only=False, turbo=True)

# plot potential
#plt.plot(xaxis, Vpot)

# plot ground state
plt.scatter(xaxis, psi, s = 1)
plt.show()


Two questions

1 - If I just put the potential (Vpot) to zero I get the classical particle in a hole solution. The solution clearly resembles the analytical one, however I get the signs mixed somehow I am not sure why this is happening.

2 - When I introduce the double level in the potential, I start getting something wrong, even here I don't understand what's happening. Any ideas?

• I’m voting to close this question because it's ultimately a code review question rather than a physics question. Dec 1, 2021 at 13:01

You could be doing the same mistake as here, On using Python to solve Time Independent Schrodinger Equation, the eigenfunctions have their values "pushed" to one of the boundaries?

With Python routines typically the columns are the eigenvectors and not the rows. That means you have to dereference your eigenvectormatrix like this [:,0] to obtain the first column/eigenvector. Derefencing by  gives you the first row and not the column.

The sign of the eigenvectors also doesn't matter. The sign/phase is just a convention but any eigenvector times a constant remains an eigenvector. The python routines do not always return the eigenvectors with the same sign. Real normed eigenvectors still have an underdetermined sign of $$\pm1$$. You can fix it by imposing a convention, for example requiring that the first coefficient is positive or picking the coefficient with highest magnitude as positive.

• Thank you very much for your reply. I get the sign issue, however I am still unsure how to fix it since it crosses zero and changes sign, but this is actually a minor concern. Regarding the column/row issue, your solution solves it partially (the points with value 0 disappear from the left half) however I still get the right half equals 0, which can't be the correct solution. Nov 30, 2021 at 11:42
• @Andrea You have flipped the sign of the whole equation and your returned eigenvalues and eigenstates are returned ordered from lowest eigenvalue to highest eigenvalues. But due to your signflip, the lowest state is at the highest energy. Try plotting [:,-1] instead of  and succesively [:,-2], [:,-3]. Those states should gives you what you would expect. I would recommend to keep the signs as in the "normal" time independent Schrödinger equation to avoid this. So that the positive E is on the right hand side. Nov 30, 2021 at 14:15
• This worked, thank you very much! Dec 1, 2021 at 10:37

First of all, I suspect one can solve this problem analytically.

As for scipy.linalg.eigh() (https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eigh.html), you just include V in matrix a.

• Thank you very much. Nov 30, 2021 at 10:59