function dYdt = kepler_ode_khan(t,y)
% assume 'Y' is a vector of form y = [ x ;  y ; vx ; vy ]
%y=[rxs;rys;vxs,vys;rxp;ryp;vxp;vyp;rxa;rya;vxa;vya];
%r = [Y(1);Y(2)];
%v = [Y(3);Y(4)];
rs = [y(1) y(2)];
rp = [y(5) y(6)];
rc = [y(9) y(10)];
rc2= [y(13) y(14)];
rxs = y(1);
rys = y(2);
vxs = y(3);
vys = y(4);
rxp = y(5);
ryp = y(6);
vxp = y(7);
vyp = y(8);
rxc = y(9);
ryc = y(10);
rxc2= y(13);
ryc2= y(14);
vxc = y(11);
vyc = y(12);
vxc2= y(15);
vyc2= y(16);
vs = [vxs vys];
vc = [vxc vyc];
vp = [vxp vyp];
vc2= [vxc2 vyc2];

% use the following three variables from the main worksheet
% must declare them 'global' here and there
global mSun mPlanet G mComet

% assumes the sun is located at r=(0 0)
Fsp = - mSun * mPlanet * G / norm(rs-rp)^3 .* (rs-rp);
Fps = -1*Fsp;

Fsc= - mSun*mComet*G/norm(rs-rc)^3.*(rs-rc);
Fcs=-1*Fsc;

as=(Fsp+Fsc)./mSun;  %acceleration of Sun
axs=as(1);
ays=as(2);

Fpc = -mPlanet*mComet*G/norm(rp-rc)^3.*(rp-rc);
Fcp=-1*Fpc;

ap=(Fps+Fpc)./mPlanet;   %acceleration of planet
axp=ap(1);
ayp=ap(2);

ac=(Fcs+Fcp)./mComet;  %acceleration of comet
axc=ac(1);
ayc=ac(2);

Fsc2= - mSun*mComet*G/norm(rs-rc2)^3.*(rs-rc2);
Fc2s=-1*Fsc2;

Fpc2 = -mPlanet*mComet*G/norm(rp-rc2)^3.*(rp-rc2);
Fc2p=-1*Fpc2;

ac2=(Fc2s+Fc2p)./mComet;  %acceleration of second comet
axc2=ac2(1);
ayc2=ac2(2);



dYdt = [vs(1); vs(2); as(1); as(2); vp(1); vp(2); ap(1); ap(2); vc(1); vc(2); ac(1); ac(2); vc2(1); vc2(2); ac2(1); ac2(2)];
end