clear
mu1 = 1;
e1=0.1;
aa1=1;
x01=aa1*(1-e1);
v01=sqrt(mu1*(1+e1)/(1-e1)/aa1);

mu2 = 0.5;
e2=0.2;
aa2=0.5;
x02=aa2*(1-e2);
v02=sqrt(mu2*(1+e2)/(1-e2)/aa2);

% Parameters for scheme
a=0;
b=100;
N=5000;  % Velocity Verlet
h=(b-a)/N;

x=zeros(N+1,2);
y=zeros(N+1,2);
u=zeros(N+1,2);
v=zeros(N+1,2);

% Initial Condition
t(1)=0;
x(1,1)=x01;
u(1,1)=0;
y(1,1)=0;
v(1,1)=v01;

x(1,2)=x02;
u(1,2)=0;
y(1,2)=0;
v(1,2)=v02;

figure(1)
axis([-1.5,1.5,-1.5,1.5])
hold on

 for i=2:N+1
    alpha1=mu1/(x(i-1,1).^2 + y(i-1,1).^2).^(3/2);
    alpha2=mu2/(x(i-1,2).^2 + y(i-1,2).^2).^(3/2);

    x(i,1)=x(i-1,1)+h*u(i-1,1)-h^2/2*alpha1*x(i-1,1);
    y(i,1)=y(i-1,1)+h*v(i-1,1)-h^2/2*alpha1*y(i-1,1);
    x(i,2)=x(i-1,2)+h*u(i-1,2)-h^2/2*alpha2*x(i-1,2);
    y(i,2)=y(i-1,2)+h*v(i-1,2)-h^2/2*alpha2*y(i-1,2);

    u(i,1)=u(i-1,1)-h/2*alpha1*x(i-1,1);
    v(i,1)=v(i-1,1)-h/2*alpha1*y(i-1,1);
    u(i,2)=u(i-1,2)-h/2*alpha2*x(i-1,2);
    v(i,2)=v(i-1,2)-h/2*alpha2*y(i-1,2);

    alpha1=1/(x(i,1).^2 + y(i,1).^2).^(3/2);
    alpha2=1/(x(i,2).^2 + y(i,2).^2).^(3/2);

    u(i,1)=u(i,1)-h/2*alpha1*x(i,1);
    v(i,1)=v(i,1)-h/2*alpha1*y(i,1);
    u(i,2)=u(i,2)-h/2*alpha2*x(i,2);
    v(i,2)=v(i,2)-h/2*alpha2*y(i,2);

    t(i)=t(i-1)+h;
    
    plot(x(i,1),y(i,1),'k.')
     plot(x(i,2),y(i,2),'b.')
    drawnow
 end
 hold off