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