function dy = orbit(x,y)

%   xspan = [0 17.0652165601579625588917206249];
%   y0 = [0.994; 0; 0; -2.00158510637908252240537862224];

mu = 0.012277471;
mustar = 1 - mu;

d1 = ((y(1) + mu)^2 + y(2)^2)^(3/2);
d2 = ((y(1)-mustar)^2+y(2)^2)^(3/2);

dy =   [ y(3)
         y(4)
         y(1)+2*y(4)-mustar*(y(1)+mu)/d1-mu*(y(1)-mustar)/d2;
         y(2)-2*y(3)-mustar*(y(2)/d1)-mu*(y(2)/d2) ];