function dz = threebodyf(t,z)
dz = zeros(4,1);
mu = 1;

r=(z(1).^2 + z(2).^2).^(3/2);
alpha=mu/r;
dz(1) = z(3);
dz(2) = z(4);
dz(3) = -alpha*z(1);
dz(4) = -alpha*z(2);

