function U0 = hillstrob(a, e)
  param = [a, e];
  U01 = lsode("hill", [1, 0, param], [0, 2*pi]);
  U02 = lsode("hill", [0, 1, param], [0, 2*pi]);
  U0 = [U01(2, 1:2)', U02(2, 1:2)'];
endfunction