function dYdt = Apollofunc(t,Y)

global mEarth mMoon mSatellite mRocket G

rm = [Y(1);Y(2);Y(3)];
vm = [Y(4);Y(5);Y(6)];

re = [Y(7);Y(8);Y(9)];
ve = [Y(10);Y(11);Y(12)];

rs = [Y(13);Y(14);Y(15)];
vs = [Y(16);Y(17);Y(18)];

rr = [Y(19); Y(20); Y(21)];
vr = [Y(22); Y(23); Y(24)];



%distance from planet to asteriod
dem = sqrt((rm(1)-re(1))^2 + (rm(2)-re(2))^2 + (rm(3)-re(3))^2);
dsm = sqrt((rs(1)-rm(1))^2 + (rs(2)-rm(2))^2 + (rs(3)-rm(3))^2);
dse = sqrt((rs(1)-re(1))^2 + (rs(2)-re(2))^2 + (rs(3)-re(3))^2);
dre = sqrt((rr(1)-re(1))^2 + (rr(2)-re(2))^2 + (rr(3)-re(3))^2);
drm = sqrt((rr(1)-rm(1))^2 + (rr(2)-rm(2))^2 + (rr(3)-rm(3))^2);

Fme = - mEarth * mMoon * G / dem^3 .* (rm-re);
Fem = - mEarth * mMoon * G / dem^3 .* (re-rm);
Fsm = - mMoon * mSatellite * G / dsm ^3 .* (rs-rm);
Fse = - mEarth * mSatellite * G / dse ^3 .* (rs-re);
Fre = - mEarth * mSatellite * G / dre ^3 .* (rr-re);
Frm = - mMoon * mSatellite * G / drm ^3 .* (rr-rm);


amx = Fme(1) / mMoon;
amy = Fme(2) / mMoon;
amz = Fme(3) / mMoon;
aex = Fem(1) / mEarth;
aey = Fem(2) / mEarth;
aez = Fem(3) / mEarth;
asx = (Fsm(1) + Fse(1))/ mSatellite;
asy = (Fsm(2) + Fse(2))/ mSatellite;
asz = (Fsm(3) + Fse(3))/ mSatellite;
arx = (Fre(1) + Frm(1))/ mRocket;
ary = (Fre(2) + Frm(2))/ mRocket;
arz = (Fre(3) + Frm(3))/ mRocket;

dYdt = [ vm(1) ; vm(2) ; vm(3); amx ; amy; amz; ve(1); ve(2); ve(3); aex; aey; aez; vs(1); vs(2); vs(3); asx; asy; asz; vr(1); vr(2); vr(3); arx; ary; arz ];  % must column vector 

