function dy=keplerf(t,y)
global gm mkr ro0 rox cs;
rr=sqrt(y(1)^2+y(2)^2);
vv=sqrt(y(3)^2+y(4)^2);
h=rr-mkr;
ro=ro0*exp(-rox*h);
dy=zeros(4,1);
dy(1)=y(3);
dy(2)=y(4);
dy(3)=-y(1)*gm/rr^3-cs*ro*vv*y(3);
dy(4)=-y(2)*gm/rr^3-cs*ro*vv*y(4);
end
