## This script illustrates the instability of the solution
## with some particular values of e and a
time = linspace(0, 50, 100);
a = .25;
e = .2;
sol = lsode("hill", [1, 0, a, e], time);


title("The solution in the phase space");
xlabel("Position");
ylabel("Velocity");
plot(sol(:,1),sol(:,2))

pause(3);
title("The position as function of time");
xlabel("Time");
ylabel("Position");
plot(time, sol(:,1));

pause(3);

clearplot;