function xdot = threebody(t,x) %Threebody ODE system %Inputs are scaler t for time and a 12x1 %vector, x, for initial positions and velocitys %of the three masses % Set values for Gravity and the three masses G = 2; m1 = 2; m2 = 2; m3 = .001; % declare and initialize 12x1 column storage xdot = zeros(12,1); x1 = x(1); u1 = x(2); x2 = x(3); u2 = x(4); x3 = x(5); u3 = x(6); y1 = x(7); v1 = x(8); y2 = x(9); v2 = x(10); y3 = x(11); v3 = x(12); %Define distances between each of the three masses r12 = sqrt( (x1-x2)^2 + (y1-y2)^2 ); r23 = sqrt( (x2-x3)^2 + (y2-y3)^2 ); r31 = sqrt( (x3-x1)^2 + (y3-y1)^2 ); %Define the ODE system xdot(1) = x(2); xdot(2) = G*m2*(x2-x1)/r12^3+G*m3*(x3-x1)/r23^3; xdot(3) = x(4); xdot(4) = G*m1*(x1-x2)/r12^3+G*m3*(x3-x2)/r23^3; xdot(5) = x(6); xdot(6) = G*m1*(x1-x3)/r31^3+G*m2*(x2-x3)/r23^3; xdot(7) = x(8); xdot(8) = G*m2*(y2-y1)/r12^3+G*m3*(y3-y1)/r31^3; xdot(9) = x(10); xdot(10) = G*m1*(y1-y2)/r12^3+G*m3*(y3-y2)/r23^3; xdot(11) = x(12); xdot(12) = G*m1*(y1-y3)/r31^3+G*m2*(y2-y3)/r23^3; % end of function