function dydt = binaryandplanet_kepler(t,y)

global ms1 ms2 mp1 G

% decompile vector y

rs1 = y(1:2);
vs1 = y(3:4);
rs2 = y(5:6);
vs2 = y(7:8);
rp1 = y(9:10);
vp1 = y(11:12);

% get distances between objects

ds1s2 = rs1-rs2;
ds1p1 = rs1-rp1;
ds2p1 = rs2-rp1;

% calculate the forces between the objects

Fs1s2 = -G*ms1*ms2/norm(ds1s2)^3 .* ds1s2;
Fs1p1 = -G*ms1*mp1/norm(ds1p1)^3 .* ds1p1;
Fs2p1 = -G*ms2*mp1/norm(ds2p1)^3 .* ds2p1;

% calculate accelerations of the objects

as1 = (Fs1s2+Fs1p1)./ms1;
as2 = (-Fs1s2-Fs2p1)./ms2;
ap1 = (-Fs1p1-Fs2p1)./mp1;

% compile dydt

dydt = [vs1(1) ; vs1(2) ; as1(1) ; as1(2) ; vs2(1) ; vs2(2) ; as2(1) ; as2(2) ; vp1(1) ; vp1(2) ; ap1(1) ; ap1(2) ];