function kepler;
global gamma mm1 mm2;
gamma=6.67*10^-11
mm1=5.97219*10^24
mm2=7.3477*10^22
r0=[0.;6400000.;0.;384400000;11060.;0.;1000.;0.];
options = odeset('RelTol',1e-7,'AbsTol',1e-7,'MaxStep',100.,'InitialStep',100.);
[t,r]=ode23(@rpp,[0 2000000], r0,options);
comet (r(:,3),r(:,4))

odeget(options,'RelTol')

function f=rpp(t,r)
global gamma mm1 mm2;
f=[r(5);r(6);r(7);r(8);
    -r(1)*gamma*mm1/sqrt(r(1)^2+r(2)^2)^3-(r(1)-r(3))*gamma*mm2/sqrt((r(1)-r(3))^2+(r(2)-r(4))^2)^3;
    -r(2)*gamma*mm1/sqrt(r(1)^2+r(2)^2)^3-(r(2)-r(4))*gamma*mm2/sqrt((r(1)-r(3))^2+(r(2)-r(4))^2)^3;
    -r(3)*gamma*mm1/sqrt(r(3)^2+r(4)^2)^3;-r(4)*gamma*mm1/sqrt(r(3)^2+r(4)^2)^3]
