function Xd = f(t, X)
    M = [1.98e30; 5.98e24; 7.35e22]; % Weights for Sun, Earth & Moon
    n = 3;           % 3 Bodies
    G = 6.67428e-11; % Gravitational constant

    Xd = zeros(4*n, 1);

    % x'(x) = v
    Xd(1:2*n) = X(2*n+1:4*n);
    
    % v'(x) = a
    for i=1:3
        for j=i+1:n
            % Compute force
            posBodyA = X(2*i-1:2*i);
            posBodyB = X(2*j-1:2*j);
            distance = norm(posBodyA - posBodyB);
            forceMagnitude = G*(M(i)*M(j)) / (distance*distance);
            forceDir = (posBodyA - posBodyB) / distance;

            % a = f/m, equal and opposite
            Xd(2*(n+i)-1:2*(n+i)) = Xd(2*(n+i)-1:2*(n+i)) ...
                   -forceDir * forceMagnitude / M(i);
            Xd(2*(n+j)-1:2*(n+j)) = Xd(2*(n+j)-1:2*(n+j)) ...
                   + forceDir * forceMagnitude / M(j);
        end
    end
end
