function ydot=anharmonic(t,y)
ydot=[y(2);-y(1)-y(1)^3];
end