function yp = dynamics2body (t, y)
  % RHS for the 1st-order system form of the 2-body problem
  yp = zeros (size (y));
  yp(1) = y(2);
  yp(2) = -y(1) / (y(1)^2 + y(3)^2)^(3/2);
  yp(3) = y(4);
  yp(4) = -y(3) / (y(1)^2 + y(3)^2)^(3/2);
endfunction