function dz = twobody2(t,z)
dz = zeros(8,1);
G = 1;
m1 = .1;
m2 = 2;
r=((z(1) - z(3)).^2 + (z(2) - z(4)).^2).^(3/2);
alpha=G/r;
dz(1) = z(5);
dz(2) = z(6);
dz(3) = z(7);
dz(4) = z(8);
dz(5) = alpha*m2*(z(3) - z(1));
dz(6) = alpha*m2*(z(4) - z(2));
dz(7) = alpha*m1*(z(1) - z(3));
dz(8) = alpha*m1*(z(2) - z(4));

% Sample entries
%[t,z] = ode45('twobody2',[0 500], [-1 0 2 0 0 -1 0 1]);
%plot(z(:,1),z(:,2),'r',z(:,3),z(:,4),'b');

% For G=1, m1=.1, m2=2 - get m1 moving around m2 and both moving 
%[t,z] = ode45('twobody2',[0 50], [-1 0 0 0 0 -1 0 0]);
%plot(z(:,1),z(:,2),'r',z(:,3),z(:,4),'b');

% phase diagram?
%  plot(z(:,1),z(:,5));