function dy=keplerf(t,y)
global g_maa r_maa ro_ohk_0 param1 param2 ;
ra=sqrt(y(1)^2+y(2)^2);
v_abs=sqrt(y(3)^2+y(4)^2);
h=ra-r_maa;
ro_ohk=ro_ohk_0*exp(-param1*h);
dy=zeros(4,1);
dy(1)=y(3);
dy(2)=y(4);
dy(3)=-y(1)*g_maa/ra^3-param2*ro_ohk*v_abs*y(3);
dy(4)=-y(2)*g_maa/ra^3-param2*ro_ohk*v_abs*y(4);

end
