Why isn't my third body plotting?

 Hi,
When I run this program, it only plots BM1 and BM2 and exlcudes BM3 and I'm completely at a loss as to why. I don't receive any errors when the program is run and the code looks alright to me. If anyone could show me as to why my third body (BM3) isn't plotting it would be greatly appreciated. Thank you.
My code is as folllows:

clc
clear
%Set Mass for Body 1
BM1 = 1*10^13;
%Set Mass for Body 2
BM2 = 10;
%Set the Mass for Body 3
BM3 = 20;
%Set the Gravitational Constant
G = 6.673*10^-11;
%Set a value fot the change in Time
DeltaT = 0.1;
%Set vales for the total time of the bodies
Total_Time = 0;
Total_Time2 = 0;
%Set the Velocity values for Body One
Vi1 = [0,0];
%Set the Velocity Values for Body Two
Vi2 = [0,7.5];
%Set the velocity values for Body Three by implementing an array
Vi3 = [0, 35];
%Set the x and y values of Si for Body 1 by implementing an array
Si1 = [0,0];
%Set the X and Y values for the initial position of Body 2 by implementing
%an array
Si2 = [10,0];
%Set the initial position (x an y values) of BM3 implementing an array
Si3 = [20,0];
%Create a while loop to calculate for the instance where ti < 10
while Total_Time < 10;
    %Calculate a radius value for use in the gravitational force equation
r1 = sqrt(((Si2(1)-Si1(1)).^2)+((Si2(2)-Si1(2)).^2));
%Calculate a radius for Bodies 3 and 2 for use in the gravitational force equation
r2 = sqrt((Si2(1)-Si2(1)).^2)+((Si3(2)-Si1(2)).^2);
%Calculate a Unit Vector to give the vector of B1 to B2 a direction
UnitV1 = (((Si1)-Si2)/(r1));
%Calculate a Unit Vector to give the vector of B1 to B2 a direction
UnitV2 = (((Si2)-Si3)/(r2)); 
%Calculate a Unit Vector for B1 to B3
%Calculate a Unit Vector for B2 to B3
%Calaculte the Gravitational force before multiplying it by the Unit Vector
F12 = (G*((BM1*BM2)/(r1.^2)));
F23 = (G*((BM2*BM3)/(r2.^2)));
%Multiply F1&2 by the Unit Vector to attain the value for F at an instant
F12 = UnitV1.*F12;
%Multiply F2&3 by the Unit Vector to attain the value for F at an instant
F23 = UnitV2.*F23;
%Calculate the velocity of the three bodies at an instant
Vi1 = -((F12.*DeltaT)./(BM1))+Vi1;
Vi2 = ((F12.*DeltaT)./(BM2))+Vi2;
Vi3 = ((F23.*DeltaT)./(BM3))+Vi3;
%Calculate the position of the two bodies at a particular instant
Si1 = (Vi1.*DeltaT)+Si1;
Si2 = (Vi2.*DeltaT)+Si2;
Si3 = (Vi3.*DeltaT)+Si3;
%Plot the Position X and Y values whil ti < 10
Total_Time = DeltaT + Total_Time;
hold all;
%Set the axis on which the orbit simulation is to be plotted
axis([-10 10 -10 10]);
%Graph the two bodies by using the scatter3 command to generate an
%isometric view
scatter3(Si1(1),Si1(2),10,'r*');
scatter3(Si2(1),Si2(2),10,'g*');
scatter3(Si3(1),Si3(2),10,'r*');
%Pause the simulation every 0.1 seconds to show the orbit of Body 2
pause(0.1)
%Use view to enable the 3D projection to be seen
view(40,35)
end

Answer 1

[Mischa Kim]

Thomas, r2 is equals zero, at least at the beginning of the loop. This causes F23 and Si3 to blow up.

It looks as if

 r2 = sqrt((Si2(1)-Si2(1)).^2)+((Si3(2)-Si1(2)).^2);

should rather be

 r2 = sqrt((Si3(1)-Si2(1)).^2)+((Si3(2)-Si1(2)).^2);

Also, you might want to adjust the axis. For starters comment out

 axis([-10 10 -10 10]);

In general, I recommend using break points for debugging.

[Thomas Stokes]

Thanks a lot for your help! However, when plotted now, the points only follow a straight line. I believe it is meant to resemble an ellipse. Would you have any recommendations as to how to achieve this elliptical shape?

[Mischa Kim]

BM2 and BM3 are very light bodies (small mass). Therefore the gravitational force between those two bodies is virtually non-existent (extremely small, that is).

 F23 = (G*((BM2*BM3)/(r2.^2)));

Hence the straight line.