function kepler_1;
global gamk;
gamma=6.67*10^-11
mkm=5.97219*10^24
mkr=6400000.
gamk=gamma*mkm/mkr^3
r0=[0.;1.1;7800/mkr;0.];
options = odeset('RelTol',1e-5,'AbsTol',1e-5,'MaxStep',1.,'InitialStep',10.);
[t,r]=ode23(@rpp,[0 10000], r0,options);
comet (r(:,1),r(:,2),0.0)
axis([-2 2. -2. 2.]);
odeget(options,'RelTol')

function f=rpp(t,r)
global gamk;
f=[r(3);r(4);-r(1)*gamk/sqrt((r(1)^2+r(2)^2))^3;-r(2)*gamk/sqrt((r(1)^2+r(2)^2))^3]
