On using Python to solve Time Independent Schrodinger Equation, the eigenfunctions have their values "pushed" to one of the boundaries?

Question

I am having trouble using numerical methods to solve Time Independent Schrodinger Equation. I am considering a quartic potential function: $$ V(x) = x^4 -4x^2.$$ $$ -\frac{d^2\psi(x)}{dx^2} + V(x) \psi(x) = E \psi(x) $$

I wish to get a few solutions of the eigenproblem (about 150).

Here is the code I have written:

import numpy as np
import matplotlib.pyplot as plt

def solveTISE(xmin,xmax,h,potFunc,params):
    xvec = np.arange(xmin,xmax,h) 
    Nx = len(xvec)
    potFunctVec = potFunc(xvec,params)
    mat = -((np.tri(Nx,Nx,1)-np.tri(Nx,Nx,-2)-3.*np.eye(Nx))/h**2-np.diag(potFunctVec))
    print(mat)
    eigenValues,eigenVectors = np.linalg.eig(mat)
    idx = eigenValues.argsort()[::1]   
    eigenValues = eigenValues[idx]
    eigenVectors = eigenVectors[:,idx]
    return eigenValues,eigenVectors,xvec
   
def doubleWell(x,params):
    return np.piecewise(x, [x < 0, x >= 0], [lambda x : params[2]*(x**4 - params[0]*x**2), lambda x :  params[3]*(x**4 - params[1]*x**2)] )

betaValue = 0.01;
hValue = 0.01;
eigEnergy1i, eigFunc1i, xvec1i = solveTISE(-10,10,hValue,doubleWell,[4,4,1,1]);

plt.plot(xvec1i,eigFunc1i[0]);
plt.show()
plt.plot(xvec1i,eigFunc1i[1]);
plt.show()

However, I get the following output:

When I changed the limits of x to : -20 and 20, the resulting wavefunction (the 0th index one/lowest energy state) looked like this:

I don't know why the wavefunction is being "pushed" to the right boundary.

Answer 1

Your script looks fine on first glance. Your problem is that your eigenfunctions are the columns of your eigFunc1 Matrix and not the rows. Use this to plot them.

plt.plot(xvec1i,eigFunc1i[:,n]);

where n is the integer of the corresponding eigenfunction. I.e. to plot the groundstate use

plt.plot(xvec1i,eigFunc1i[:,0]);

Also note that some_array[::1] doesn't change the array, as you you are stepping over all values with a step size of 1.

EDIT:

Due to popular demand a small plot of the first eigenvectors n=0 and n=1 at the height of their eigenvalue as baseline which can be read off at the left hand y-axis.

The figure was created with

fig, ax = plt.subplots()
ax.plot(xvec1i,eigFunc1i[:,0]+eigEnergy1i[0], c="b", lw=3);
ax.plot(xvec1i,eigFunc1i[:,1]+eigEnergy1i[1], c="orange", lw=3);
ax.axhline(eigEnergy1i[0],c="b", alpha=0.5, lw=3, label="n=0")
ax.axhline(eigEnergy1i[1], c="orange", alpha=.5, lw=3, label="n=1")
ax_twinx = ax.twinx()
ax_twinx.plot(xvec1i,doubleWell(xvec1i,[4,4,1,1]))
ax.legend()
ax.set_xlim(-10, 10)

Answer 2

Your another (small) problem is that your values of $x$, called xvec, are not in the range $\{-10, -9.99, \dots, +9.99, +10\}$. Instead they lack the $+10$, because arange(begin,end,step) returns the range [begin,end). The result for your code is that the states that should have very close energies will instead differ more than expected, due to the broken symmetry. The fix for this is to replace your generation of xvec with

    xvec = np.arange(xmin,xmax+h,h) 

This is much better, but it still gives you not ideal values, so the precision of the eigenvalues will be a bit smaller than you may like. Here's the output of print(xvec) after setting np.set_printoptions(precision=17):

[-10.                 -9.99               -9.98              ...
   9.979999999999574   9.989999999999574   9.999999999999574]

I suppose it's due to the way arange() is implemented. To fix this, you can try to make sure that the step is integral, so that you don't accumulate too much rounding error. An example (just a dirty hack for demonstration, not intended for production):

    xvec = np.arange(xmin*1000,(xmax+h)*1000,h*1000)/1000

This code does generate

[-10.    -9.99  -9.98 ...   9.98   9.99  10.  ]