Python: odeint to solve ODEs with variable coefficients (QHO)

Question

I'm trying to solve the equation y'' + (epsilon-x^2)y = 0 numerically using odeint. I know the solutions (the wavefunctions of a QHO), but the output from odeint has no apparent relation to it. I can solve ODEs with constant coefficients just fine, but as soon as I move to variable ones, I can't solve any of the ones I tried. Here's my code:

#!/usr/bin/python2
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as spi

x = np.linspace(-5,5,1e4)

n = 0
epsilon = 2*n+1 

def D(Y,x):
    return np.array([Y[1], (epsilon-x**2)*Y[0]])

Y0 = [0,1]

Y = spi.odeint(D,Y0,x)
# Y is an array with the first column being y(x) and the second y'(x) for all x
plt.plot(x,Y[:,0],label='num')
#plt.plot(x,Y[:,1],label='numderiv')

plt.legend()
plt.show()

And the plot: [not enough rep:] https://drive.google.com/file/d/0B6840LH2NhNpdUVucUxzUGFpZUk/edit?usp=sharing

Look here for plots of solution: http://hyperphysics.phy-astr.gsu.edu/hbase/quantum/hosc5.html

Answer 1

It looks like your equation is not correctly interpreted. You have a differential equation y'' + (epsilon-x^2)y = 0, but you forget a minus sign in your vector form. In particular it should be

y[0]' = y[1]
y[1]' = -(epsilon-x^2)y[0]

So (adding the minus sign in front of the epsilon term

def D(Y,x):
    return np.array([Y[1], -(epsilon-x**2)*Y[0]])

In fact the plot you have is consistent with the DE y'' + (epsilon-x^2)y = 0. Check it out: Wolphram Alpha