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.