function dYdt = kepler_ode_khan2(t,y_new)
% 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_new(1) y_new(2)];
rp = [y_new(5) y_new(6)];
rc3 = [y_new(9) y_new(10)];
rc4= [y_new(13) y_new(14)];
rxs = y_new(1);
rys = y_new(2);
vxs = y_new(3);
vys = y_new(4);
rxp = y_new(5);
ryp = y_new(6);
vxp = y_new(7);
vyp = y_new(8);
rxc3 = y_new(9);
ryc3 = y_new(10);
rxc4= y_new(13);
ryc4= y_new(14);
vxc3 = y_new(11);
vyc3 = y_new(12);
vxc4= y_new(15);
vyc4= y_new(16);
vs = [vxs vys];
vc3 = [vxc3 vyc3];
vp = [vxp vyp];
vc4= [vxc4 vyc4];

% 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;

Fsc3= - mSun*mComet*G/norm(rs-rc3)^3.*(rs-rc3);
Fc3s=-1*Fsc3;

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

Fpc3 = -mPlanet*mComet*G/norm(rp-rc3)^3.*(rp-rc3);
Fc3p=-1*Fpc3;

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

ac3=(Fc3s+Fc3p)./mComet;  %acceleration of comet
axc3=ac3(1);
ayc3=ac3(2);

Fsc4= - mSun*mComet*G/norm(rs-rc4)^3.*(rs-rc4);
Fc4s=-1*Fsc4;

Fpc4 = -mPlanet*mComet*G/norm(rp-rc4)^3.*(rp-rc4);
Fc4p=-1*Fpc4;

ac4=(Fc4s+Fc4p)./mComet;  %acceleration of second comet
axc4=ac4(1);
ayc4=ac4(2);



dYdt = [vs(1); vs(2); as(1); as(2); vp(1); vp(2); ap(1); ap(2); vc3(1); vc3(2); ac3(1); ac3(2); vc4(1); vc4(2); ac4(1); ac4(2)];
end