function yp=Arenstorf(t,y);
% ARENSTORF ODE right hand side for the Arenstorf orbit problem
%   yp=Arenstorf(t,y); describes a system of ODEs which model the
%   flight of a light object between the earth and the moon.

a=0.012277471; b=1-a;
D1=((y(1)+a)^2+y(2)^2)^(3/2);
D2=((y(1)-b)^2+y(2)^2)^(3/2);
yp(1,1)=y(3);
yp(2,1)=y(4);
yp(3,1)=y(1)+2*y(4)-b*(y(1)+a)/D1-a*(y(1)-b)/D2;
yp(4,1)=y(2)-2*y(3)-b*y(2)/D1-a*y(2)/D2;
yp=yp(:);                                 % make a column vector