function dw=threebody(t,w)  
%Computes changes in y for three body problem.
%  The vector w is partioned as:
%  w(1,2)=(x1,y1)  Position of body 1
%  w(3,4)=(x2,y2)  Position of body 2
%  w(5,6)=(x3,y3)  Position of body 3
%  w(7,8)=(dx1,dy1) Velocity of body 1
%  w(9,10)=(dx2,dy2) Velocity of body 2
%  w(11,12)=(dx3,dy3) Velocity of body 3

%G is gravity, m1, m2, m3 are the masses...
G=1;
m1=5;m2=1;m3=1;

%The softening constant ensures no blowups!
soft=0.2;

dw(1:6)=w(7:12);

%Distance computations
r12=(norm(w(1:2)-w(3:4))^2+soft^2)^1.5;
r13=(norm(w(1:2)-w(5:6))^2+soft^2)^1.5;
r23=(norm(w(3:4)-w(5:6))^2+soft^2)^1.5;

dw(7:8)=((G*m2)/r12)*(w(3:4)-w(1:2))+((G*m3)/r13)*(w(5:6)-w(1:2));
dw(9:10)=((G*m1)/r12)*(w(1:2)-w(3:4))+((G*m3)/r23)*(w(5:6)-w(3:4));
dw(11:12)=((G*m1)/r13)*(w(1:2)-w(5:6))+((G*m3)/r23)*(w(3:4)-w(5:6));
dw=dw(:);