  function du = Kepler(t,u)
% Odefunction for 2-body problem (C. van Loan p. 343)
% du = Kepler(t,u)  % 
% t (time) is a scalar and u is a 4-vector whose components satisfy
%
%            u(1) = x      u(2) = x'
%            u(3) = y      u(4) = y'
%
% where (x(t),y(t)) is the position of the satellite w.r.t.Earth
% du is a 4-vector that is the derivative of u at time t.

r3 = (u(1)^2 + u(3)^2)^(3/2);
du = [  u(2)     ;
       -u(1)/r3  ;
        u(4)     ;
       -u(3)/r3] ;