# function kepler useful for solving
# planetary motion around sun.

function xdot = kepler(x,t)
  xdot = zeros(4,1);
# GM_sun = 4*pi^2
# when [time] = yr and [r] = AU
a = 39.5;
xdot(1) = x(2);
xdot(2) = -a*x(1)/((x(1))**2 + (x(3))**2)^1.5;
xdot(3) = x(4);
xdot(4) = -a*x(3)/((x(1))**2 + (x(3))**2)^1.5;

endfunction

