function forcee=F_body(X)
        global G
        forcee=zeros(3,3,3);
        for i=1:3
            for j=1:3
                for k=1:3
                    distances=dis(X);
                    if i~=j
                forcee(i,j,k)=(G*X(7,i)*X(7,j)/distances(i,j)^3)*(X(k,j)-X(k,i));
                    else if i==j
                            forcee(i,j,k)=0;
                        end
                    end
                end
            end
        end

Published with MATLAB® R2019a