function Y = grav2(t,x,m)

Y = zeros(12,1);
% coupled differential equations
% x = [x,y,vx,vy]

% 3D orbits
for i=1:2
    % dx/dt = v
    Y(6*i-3, 1) = x(6*i);
    Y(6*i-4, 1) = x(6*i-1);
    Y(6*i-5, 1) = x(6*i-2);
    
    a = 0;
    b = 0;
    c = 0;
    for j=1:2
        if i~=j
            %need to work out the correct order and subscripts
            r = sqrt((x(6*j-3)-x(6*i-3))^2+(x(6*j-4)-x(6*i-4))^2+(x(6*j-5)-x(6*i-5))^2);
            
            % dv/dt = a  (a = F/m) (F = GMm|x2-x1|/r^3)
            
            a = m(i)*(x(6*i-5)-x(6*j-5))/r^3;
            b = m(i)*(x(6*i-4)-x(6*j-4))/r^3;
            c = m(i)*(x(6*i-3)-x(6*j-3))/r^3;
            
            Y(6*j-2, 1) = Y(6*j-2, 1) + a;
            Y(6*j-1, 1) = Y(6*j-1, 1) + b;
            Y(6*j, 1) = Y(6*j, 1) + c;
        end
    end
end