global gamk;
gamma=6.67402*10^-11;
m_maa=5.97219*10^24;
mkr=6378.;
gamk=gamma*m_maa*10^-9;
r0=[0.;6378.;10;0.];
options = odeset('RelTol',1e-5,'AbsTol',1e-5,'MaxStep',1.,'InitialStep',1.);
[t,r]=ode45(@kplf,[0,400000], r0,options);

hold on;
axis equal;
grid on;
z0=zeros(1,length(t));
im= imread('earth1.jpg');
[x,y,z]=sphere(50);
surf(x,y,-z,im,'FaceColor','texturemap','Edgecolor','none');

plot3(r(:,1)/mkr,r(:,2)/mkr,z0);
comet3(r(:,1)/mkr,r(:,2)/mkr,z0,0.01);

function  dy=kplf(t,r)
global gamk;
dy=zeros(4,1);
dy(1)=r(3);
dy(2)=r(4);
dy(3)=-r(1)*gamk/sqrt((r(1)^2+r(2)^2))^3;
dy(4)=-r(2)*gamk/sqrt((r(1)^2+r(2)^2))^3;
end