T=linspace(0, 50, 1000);
y = lsode("twobody", [1;1;1;1], T);
#plot(y(:,1), y(:,2))
plot(T, y(:,1))
