Code to model a spacecraft is given by two 2nd-order equations for the position of the body (x(t), y(t)) --- Numerical Method: Störmer-Verlet (Dual) -- Matlab

DualVeerlet.m

%---------------- Setting Problem --------------------%
x0  =  0.994;                % x(0)  Initial position
y0  =  0.0;                  % y(0)  Initial position
dx0 =  0.0;                  % x'(0) Initial velocity
dy0 = -2.031732629557;       % y'(0) Initial velocity
%-----------------------------------------------------%
%---------------- Setting Model ----------------------%
mu1 = 0.012277471;          % Moon mass,  M_m/(M_e+M_m) 
mu2 = 1-mu1;                % Earth mass, M_e/(M_e+M_m) 
%-----------------------------------------------------%
  D1  = @(x,y) ( (x+mu1).^2+y.^2 ).^(1.5);
  D2  = @(x,y) ( (x-mu2).^2+y.^2 ).^(1.5);
   r  = @(x,y) 1-( mu1./D2(x,y)+mu2./D1(x,y) );
   g  = @(x,y) mu1*mu2*(1./D2(x,y)-1./D1(x,y));
   f  = @(q,v) [ r(q(1),q(2))*q(1)+g(q(1),q(2))+2*v(2);...
                 r(q(1),q(2))*q(2)-2*v(1)]; 
             
   l  = @(q) [ r(q(1),q(2))*q(1)+g(q(1),q(2));...
               r(q(1),q(2))*q(2)];        
%-----------------------------------------------------%
dt = 0.0001;
t  = 0;
T  = 11.12455; % period  
N = ceil( T/dt)+2;
A = (1/(dt^2+1))*[1 dt;-dt 1];
q  = zeros(2,N); v  = zeros(2,N);
q(:,1)=[x0;y0]; v(:,1)=[dx0;dy0];
for n=2:N
    vn12 = v(:,n-1);
    qn12 = q(:,n-1);
    qn   = qn12+(dt/2)*vn12;
    v(:,n) = A*(vn12+(dt/2)*(f(qn,vn12)+l(qn)));
    q(:,n) = qn+(dt/2)*v(:,n);
end
plot(q(2,:),q(1,:),'-.b','LineWidth',2) % plot(y,x)
hold on
plot(0,-mu1,'ro',0,mu2,'bo','MarkerSize', 12)