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

