function ydot=orbit2df(t ,y)
global N G M

% set the velocities
ydot=zeros(4*N,1);
for i=1:(2*N)
    ydot(i)=y(2*N+i );
end

a1=0; 
a2=0; 
for i=1:N 
    for j=1:N
        r=(y(2*i-1)-y(2*j-1))^2+(y(2*i)-y(2*j ))^2;
        if j~=i 
            a1=a1+(M(1,j)/(sqrt(r))^3)*(y(2*j-1)-y(2*i-1));
            a2=a2+(M(1,j)/(sqrt(r))^3)*(y(2*j)-y(2*i ));
        end
    end
    ydot(2*N + 2*i-1)=G*a1;
    ydot(2*N + 2*i)=G*a2;
    a1=0;
    a2=0;
end
