function wdot = hill(w, t) 
  a = w(3);
  e = w(4);
  wdot(1) = w(2);
  wdot(2) = - (a + e * cos(t)) / (1 + e * cos(t)) * w(1);
  wdot(3) = 0;
  wdot(4) = 0;
endfunction