Simulating gravity, forces seem to be miscalculated (need help)

I am trying to simulate how 5-15 different bodies affect eachother in terms of gravity.

This is my first time coding, so needless to say I need some help.

The problem is that the forces acting on the bodies seem to be miscalculated or at least are acting in the wrong direction, for instance all forces are positive in the x-axis. Some of the bodies seem to "run away" when another gets close.

Here is my code so far (sorry if it is a mess).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%GRAVITATION %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
clc;
n=randi([5,15]);                        %n is the amount of bodies
XYM=rand(3,n)*(1000-10);                %A matrix with x- and y-coordinates and mass of each body. (Each column is a 
                                        body)
t=0;
v0X=zeros(1,n);
v0Y=zeros(1,n);
G=6.673e-11;                            

while n>1

    sz=XYM(3,:);                            
    c=gradient(XYM(3,:));                   

    scatter(XYM(1,:),XYM(2,:),sz,c,'filled')
    whitebg('k')

    pause(.02)

    %Distances are calculated, also in x and y

    m=n;
    dist=zeros(m,n);
    for i1=1:m
        for j1=1:n

            dist(j1, i1) = sqrt((XYM(1,i1) - XYM(1, j1)) ^ 2 + ...
                        (XYM(2,i1) - XYM(2,j1)) ^ 2);
        end
    end

    distX=zeros(m,n);
    for i2=1:m
        for j2=1:n

            distX(j2, i2) = XYM(1, i2) - XYM(1,j2);

        end
    end

    distY=zeros(m,n);
    for i3=1:m
        for j3=1:n

            distY(j3, i3) = XYM(2, i3) - XYM(2,j3);

        end
    end

    %The forces are calculated and summed

    F=zeros(m,n);
    for i4=1:m
        for j4=1:n

            F(j4,i4)=(G*XYM(3,i4)*XYM(3,j4))/(dist(i4,j4).^2);

        end
    end

    F(~isfinite(F))=0;

    k=distY./distX;
    k(isnan(k))=0;
    angle=atan(k);

    FX=F.*cos(angle);
    FY=F.*sin(angle);
    FX=sum(FX);
    FY=sum(FY);

    %Force (FX,FY) -> acceleration (aX,aY) -> velocity (v0X,v0Y) -> distance 
    travelled (sX,sY)

    aX=FX./XYM(3,:);
    aY=FY./XYM(3,:);

    v0X=v0X+aX*t
    v0Y=v0Y+aY*t;

    sX=v0X*t+(aX.*t^2)/2;
    sY=v0Y*t+(aY.*t^2)/2;

    XYM(1,:)=XYM(1,:)+sX;
    XYM(2,:)=XYM(2,:)+sY;

    %Bodies put together when less or equal to 10 units away from eachother

    d=dist<=10;
    d=d-eye(m,n);
    c=sum(sum(d));

    if c>0

        [r,c]=find(d);
        a=r(1,1);
        b=r(2,1);

        v0X(:,a)=(XYM(3,a).*v0X(:,a)+XYM(3,b).*v0X(:,b))/(XYM(3,a)+XYM(3,b));
        v0Y(:,a)=(XYM(3,a).*v0Y(:,a)+XYM(3,b).*v0Y(:,b))/(XYM(3,a)+XYM(3,b));

        v0X(:,b)=[];
        v0Y(:,b)=[];

        XYM(3,a)=XYM(3,a)+XYM(3,b);
        XYM(:,b)=[];

        n=n-1;
        m=n;

    end

    t=t+60;

end

What are the flaws of the code, if any? Am I even on the right track?

[Guillaume]

I also don't really have time to go through the code and try to understand it (it needs a lot more comments, each variable should say what its purpose is). You can certainly simplify your code and probably get rid of all the loops. For example, to calculate dist you can use pdist2 or simply:

%no need to preallocate dist, just calculate it all in one go with:
dist = hypot(XYM(1, :) - XYM(1, :)', XYM(2, :) - XYM(2, :)');

[Viktor Bodin]

Okay, understandable, I'll work on the comments! Thanks for the advice :)

[Viktor Bodin]

This is the updated code with comments

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%GRAVITATION %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
clc;
n=randi([5,15]);                        %n is the amount of bodies
XYM=rand(3,n)*(1000-10);                %A matrix with x- and y- 
                                        %coordinates 
                                        %and mass of each body. 
                                        %(Each column is a body)
t=0;
v0X=zeros(1,n);
v0Y=zeros(1,n);
G=6.673e-11;                            

while n>1

    sz=XYM(3,:);        
    %I tried to make the colors of each body depend on the mass of it, but
    %it doesn't seem to work
    c=gradient(XYM(3,:));                   

    scatter(XYM(1,:),XYM(2,:),sz,c,'filled')
    whitebg('k')

    pause(.02)

    %%Distances are calculated, also in x and y

    %These are the absolute distances between all bodies
    dist = hypot(XYM(1, :) - XYM(1, :)', XYM(2, :) - XYM(2, :)');       

    %These are the distances in x between all bodies
    m=n;
    distX=zeros(m,n);                                                   
    for i2=1:m
        for j2=1:n

            distX(j2, i2) = XYM(1, i2) - XYM(1,j2);

        end
    end

    %These are the distances in y between all bodies
    distY=zeros(m,n);                                                   
    for i3=1:m
        for j3=1:n

            distY(j3, i3) = XYM(2, i3) - XYM(2,j3);

        end
    end

    %Forces are calculated, also in x and y, which is why the distances in
    %x and y were needed

    F=zeros(m,n);
    for i4=1:m
        for j4=1:n

            F(j4,i4)=(G*XYM(3,i4)*XYM(3,j4))/(dist(i4,j4).^2);

        end
    end

    F(~isfinite(F))=0;

    %k is a variable used for simpler calculation of the angles
    k=distY./distX;
    k(isnan(k))=0;
    %The angles from the x-axis to the line between each body to all other 
    %bodies are calculated, which will be needed for calculating the 
    %forces in x and y 
    angle=atan(k);

    %The forces in each axis are calculated from the angles and then summed
    FX=F.*cos(angle);
    FY=F.*sin(angle);
    FX=sum(FX);
    FY=sum(FY);

    %Acceleration in x and y, which will be used to calculate 
    %velocity (a=F/m)
    aX=FX./XYM(3,:);                
    aY=FY./XYM(3,:);

    %Velocities in x and y, which will be used to calculate distance 
    %travelled (v=v0+a*t)
    v0X=v0X+aX*t;                    
    v0Y=v0Y+aY*t;

    %The distances each body travels in t amount of time in x 
    %and y (lenght=v0*t+(a*t^2)/2)
    sX=v0X*t+(aX.*t^2)/2;           
    sY=v0Y*t+(aY.*t^2)/2;

    %The distances sX and sY are added to each body's current position
    XYM(1,:)=XYM(1,:)+sX;           
    XYM(2,:)=XYM(2,:)+sY;

    %%Bodies will be put together when less or equal to 10 units away 
    %from each other

    %d is a matrix with the value 1 if bodies are less or equal to 10 
    %units from one another
    d=dist<=10;                     
    d=d-eye(m,n);

    %c sums the matrix d
    c=sum(sum(d));                  

    %if c>0 two bodies are close and need to be added
    if c>0                          

        %r will be used to find which bodies are close to one another by
        %retrieving the coordinates of each body in the XYM-matrix, with
        %the use of d
        [r,c]=find(d);
        %a is the first body in the collision, given by the coordinates 
        %provided by r
        a=r(1,1);  
        %b is the second body in the collision
        b=r(2,1);                   

        %The new velocities (in x and y) from the collision are calculated 
        %and replaces the velocities of a
        v0X(:,a)=(XYM(3,a).*v0X(:,a)+XYM(3,b).*v0X(:,b))/(XYM(3,a)+XYM(3,b));       
        v0Y(:,a)=(XYM(3,a).*v0Y(:,a)+XYM(3,b).*v0Y(:,b))/(XYM(3,a)+XYM(3,b));

        %The velocities of b are deleted
        v0X(:,b)=[];                
        v0Y(:,b)=[];

        %the masses of a and b are added and replaces the original 
        %mass of a
        XYM(3,a)=XYM(3,a)+XYM(3,b);
        %Body b is deleted
        XYM(:,b)=[];                    

        %n is the amount of bodies, now reduced
        n=n-1; 
        %m is used in loops and need to be equal to n
        m=n;                            

    end

    %60 seconds pass for each loop
    t=t+60;                             

end

Answer 1

[James Tursa]

I don't have the time to step through your code, but I don't see where you set the sign of F as an attracting force. Maybe all you need to do is flip the sign of F just before you start calculating the accelerations:

F = -F;

[Viktor Bodin]

Thanks for the reply! This didn't help too much unfortunatly :( I guess it's slightly better but still far from perfect. I might have messed up the calculations from the very beginning, not sure though.