function dYdt = kepler_n_ode(t,Y)
global G b
n = size(Y,1)/8;
X = zeros(size(Y));
for i = 0:n-1
    ii = i * 8;
    r = [Y(ii+1);Y(ii+2);Y(ii+3)];          % unpack radius vector
    v = [Y(ii+4);Y(ii+5);Y(ii+6)];          % unpack velocity vector
    m = Y(ii+7);                            % unpack mass scaler
    p = Y(ii+8);                            % unpack density scaler
    A = pi*((3*m)/(4*pi*p))^(2/3);          % calculate reference area
    a = (b*norm(v)*A/m).*(-v);              % calculate acceleration due to drag
    for j = 1:n-1
        jj = j * 8;
        if j < ( i + 1 )
            jj = ( ( j - 1 ) * 8 );
        end
        r  = [Y(jj+1) ; Y(jj+2) ; Y(jj+3)] - r;
        mm = m * Y(jj+7);                   % determine the product of the two masses
        a  = a + (G*mm/(m*norm(r)^3)).*r;   % calculate acceleration due to gravity
    end
    X(ii+1:ii+6) = [ v ; a ];
end
dYdt = X;