function xdot = twobody(x, t)
  fric = .000;
  xdot(1) = x(3);
  xdot(2) = x(4);
  rv = x(1:2);
  d = norm(rv).^3; 
  ## Less efficient alternative
  ##   d = sqrt(x(1)^2 + x(2)^2)^3
  xdot(3) = - x(1) / d - fric * x(3);
  xdot(4) = - x(2) / d - fric * x(4);
endfunction
