Problem with orbital simulation.

Hi, When I go to run this script, instead of resembling an elliptical orbit the BM3 component generates a large straight line. I was wondering why this line is generated instead of a ellipse and if there are any ways to remedy this problem. I believe it may involve vectors. My code is as follows:

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, 15];
%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 = [15,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((Si3(1)-Si2(1)).^2)+((Si3(2)-Si2(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([-100 100 -100 100]);
%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,'b*');
%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

Any help would be greatly appreciated!

(I have also attached an image of the problem)

[Joseph Cheng]

Well without digging into too much, check your masses that you specifically want to use and velocities. For instance i see that based on your relationships you are trying to have BM2 referenced to BM1. but then have a larger mass BM3 referenced to BM2 and not BM1.

If you display F23 you can see that based on the values that you've used there is little influence on the two bodies. Without knowing what you're trying to model you may want to include the force between BM1 and BM3.

[Thomas Stokes]

Thank you for your reply Joseph, but how would I include the force between BM1 and BM3?