function ydot = kepler_ode(t,y)
%obtain ?r? and ?v? from ?y? 
rs=[y(1) y(2)];
vs=[y(3) y(4)];

rp=[y(5) y(6)];
vp=[y(7) y(8)];

global ms mp G

Fs= -(G*(ms*mp)./(norm(rs-rp))^3).*(rs-rp);
%Fdrag = -(10^19)*vp; % (c=2 here is arbitrary)
Fdrag= (-(10^28)/norm(rs-rp)).*vp;

Fp=-Fs+Fdrag;

Fdrag/Fs

as = Fs/ms;
ap= Fp/mp;

ydot =[ vs  as vp  ap]';
