## Implementation in Python

Let's turn our attention to actually implementing these ideas on the computer. Here is the code [github:potentialwell.py]:

```
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,4,1000)
dx = x[1]-x[0]
kinetic_energy = -0.5/dx**2 \
* (np.diagflat(-2*np.ones(1000))\
+np.diagflat(np.ones(999),-1)\
+np.diagflat(np.ones(999),1))
energies,vec = np.linalg.eigh(kinetic_energy)
print(energies[0:4])
plt.plot(x,vec[:,0])
plt.plot(x,vec[:,1])
plt.plot(x,vec[:,2])
plt.show()
```

It really is quite simple. First, we import the necessary libraries. Then we create the x-axis (the command np.linspace(0,4,1000) means: create 1000 evenly spaced points between 0 and 4). We then make the kinetic energy matrix. Diagflat creates a matrix with the given values on the diagonal; np.ones(1000) creates an array of length 1000 filled with ones. The second argument in diagflat, -1 or 1, means "place it one level below (-1) or one level above (+1) the diagonal".

The difficult part of the code --- actually solving the eigenvalue problem --- is provided by numpy: we will not write our own eigenvalue algorithms in this book. In general, it is a waste of time to do that (except for learning purposes), because many highly optimized libraries exist for this purpose. In numpy, np.eigh(matrix) solves the eigenvalues and eigenvectors of the given (Hermitian) matrix.

Now, if you run this code, it should print out the energies as about [0.307, 1.229,2.765,4.915]. As I indicated previously, the analytical solution of this problem is known; the eigenvalues are given by

\begin{align} E_n = \frac{n^2\pi ^2}{2L^2}, \end{align}where \(L\) is the energy level and \(n\) the number of the energy level. A quick check with a calculator for \(L=4\) and \(n\in \{ 1,2,3,4\}\) confirms that our values are correct.

If we now wanted to add a potential --- say the harmonic potential \(x^2\) --- the process would be rather simple. You'd just add x**2 to the diagonal, like so:

```
hamilton = kinetic_energy + np.diagflat(x**2)
```

and then replace the matrix in np.eigh with "hamilton". So changing the potential really is quite easy.

Let's see what happens if we do set periodic boundary conditions (without the potential; if you want periodic boundary conditions, the potential should also be periodic, unlike \(x^2\)) as discussed in the previous section. The only thing that should be changed is the upper right and lower left corner of the matrix as discussed. This is easy to do; the modified code [github:periodic.py] reads

```
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,4,1000)
dx = x[1]-x[0]
kinetic_energy = -0.5/dx**2 \
* (np.diagflat(-2*np.ones(1000))\
+np.diagflat(np.ones(999),-1)\
+np.diagflat(np.ones(999),1))
kinetic_energy[0,999] = -0.5/dx**2
kinetic_energy[999,0] = -0.5/dx**2
energies,vec = np.linalg.eigh(kinetic_energy)
print(energies[0:4])
plt.plot(x,vec[:,0])
plt.plot(x,vec[:,1])
plt.plot(x,vec[:,2])
plt.show()
```

You should get different energy eigenvalues and different wave functions. You should also note that our wave function indeed is the same at the beginning and the end of he grid. Hence, we can modify the behavior at the boundaries by reasoning about the matrix elements.

It is very important to understand what is going on in this chapter, as this is the foundation on which the rest of the book is built. Hence, we now turn our attention to one final example, the Kronig-Penney model.