threebody.m

The programme ahout scientific computation

function ydot = threebody(t,y);
%   THREEBODY  Three body gravitational attraction.
%   Use with ODE solvers.
p1 = y(1:3);
p2 = y(4:6);
p3 = y(7:9);
r12 = p1-p2;
r23 = p2-p3;
r31 = p3-p1;
q1 = -r12/norm(r12)^3+r31/norm(r31)^3;
q2 = -r23/norm(r23)^3+r12/norm(r12)^3;
q3 = -r31/norm(r31)^3+r23/norm(r23)^3;
ydot = [y(10:18); q1; q2; q3];