function dYdt = kepler_ode(t,Y)

rA = [Y(1);Y(2);Y(3)];
vA = [Y(4);Y(5);Y(6)];
rS = [Y(7);Y(8);Y(9)];
vS = [Y(10);Y(11);Y(12)];
rMe = [Y(13);Y(14);Y(15)];
vMe = [Y(16);Y(17);Y(18)];
rV = [Y(19);Y(20);Y(21)];
vV = [Y(22);Y(23);Y(24)];
rE = [Y(25);Y(26);Y(27)];
vE = [Y(28);Y(29);Y(30)];
rM= [Y(31);Y(32);Y(33)];
vM= [Y(34);Y(35);Y(36)];
rJ = [Y(37);Y(38);Y(39)];
vJ = [Y(40);Y(41);Y(42)];

global mS mJ mE mA G mM mMe mV
mJ  = 10e+23+mJ;
mMe = 10e+22+mJ;
mV  = 10e+22+mJ;
mE  = 10e+22+mJ;
mM  = 10e+21+mJ;

X   = norm(rS-rA);% SUN ASTEROID
X1  = norm(rS-rMe);% SUN MERCURY
X2  = norm(rS-rV);% SUN VENUS
X3  = norm(rS-rE);% SUN EARTH
X4  = norm(rS-rM);% SUN MARS
X5  = norm(rS-rJ);% SUN JUPITER
X6  = norm(rA-rMe);% ASTEROID MERCURY
X7  = norm(rA-rV);% ASTEROID VENUS
X8  = norm(rA-rE);% ASTEROID EARTH
X9  = norm(rA-rM);% ASTEROID MARS
X10 = norm(rA-rJ);% ASTEROID JUPITER
X11 = norm(rMe-rV);% MERCURY VENUS
X12 = norm(rMe-rE);% MERCURY EARTH
X13 = norm(rMe-rM);% MERCURY MARS
X14 = norm(rMe-rJ);% MERCURY JUPITER
X15 = norm(rV-rE);% VENUS EARTH
X16 = norm(rV-rM);% VENUS MARS
X17 = norm(rV-rJ);% VENUS JUPITER
X18 = norm(rE-rM);% EARTH MARS
X19 = norm(rE-rJ);% EARTH JUPITER
X20 = norm(rM-rJ);% MARS JUPITER


% Force and Accerlation Calculations:

FA= (((-G*mA*mS)/(X^3)).*(rA-rS))+(((-G*mA*mMe)/(X6^3)).*(rA-rMe))+(((-G*mA*mV)/(X7^3)).*(rA-rV))+(((-G*mA*mE)/(X8^3)).*(rA-rE))+(((-G*mA*mM)/(X9^3)).*(rA-rM))+(((-G*mA*mJ)/((X10)^3)).*(rA-rJ));
aA= FA./mA;

FS = (((-G*mS*mA)/((X)^3)).*(rS-rA))+(((-G*mS*mMe)/((X1)^3)).*(rS-rMe))+(((-G*mS*mV)/(X2^3)).*(rS-rV))+(((-G*mS*mE)/(X3^3)).*(rS-rE))+(((-G*mS*mM)/(X4^3)).*(rS-rM))+(((-G*mS*mJ)/(X5^3)).*(rS-rJ));
aS = FS ./ mS;

FMe=(((-G*mMe*mA)/(X6^3)).*(rMe-rA))+(((-G*mMe*mS)/(X1^3)).*(rMe-rS))+(((-G*mMe*mV)/(X11^3)).*(rMe-rV))+(((-G*mMe*mE)/(X12^3)).*(rMe-rE))+(((-G*mMe*mM)/(X13^3)).*(rMe-rM))+(((-G*mMe*mJ)/(X14^3)).*(rMe-rJ));
aMe= FMe./mMe;

FV =(((-G*mV*mA)/(X7^3)).*(rV-rA))+(((-G*mV*mS)/(X2^3)).*(rV-rS))+(((-G*mV*mMe)/(X11^3)).*(rV-rMe))+(((-G*mV*mE)/(X15^3)).*(rV-rE))+(((-G*mV*mM)/(X16^3)).*(rV-rM))+(((-G*mV*mJ)/(X17^3)).*(rV-rJ));
aV = FV./mV;

FE= (((-G*mE*mA)/(X8^3)).*(rE-rA))+(((-G*mE*mS)/(X3^3)).*(rE-rS))+(((-G*mE*mMe)/(X12^3)).*(rE-rMe))+(((-G*mE*mV)/(X15^3)).*(rE-rV))+(((-G*mE*mM)/(X18^3)).*(rE-rM))+(((-G*mE*mJ)/(X19^3)).*(rE-rJ));
aE=FE./mE;

FM =(((-G*mM*mA)/(X9^3)).*(rM-rA))+(((-G*mM*mS)/(X4^3)).*(rM-rS))+(((-G*mM*mMe)/(X13^3)).*(rM-rMe))+(((-G*mM*mV)/(X16^3)).*(rM-rV))+(((-G*mM*mE)/(X18^3)).*(rM-rE))+(((-G*mM*mJ)/(X20^3)).*(rM-rJ));
aM = FM./mM;

FJ= (((-G*mJ*mA)/(X10^3)).*(rJ-rA))+(((-G*mJ*mS)/((X5)^3)).*(rJ-rS))+(((-G*mJ*mMe)/((X14)^3)).*(rJ-rMe))+(((-G*mJ*mV)/(X17^3)).*(rJ-rV))+((((-G*mJ*mE)/(X19^3)).*(rJ-rE)))+(((-G*mJ*mM)/(X20^3)).*(rJ-rM));
aJ= FJ ./mJ;
    
dYdt = [ vA(1) ; vA(2) ; vA(3); aA(1) ; aA(2) ; aA(3) ; vS(1) ; vS(2) ; vS(3) ; aS(1) ; aS(2) ; aS(3) ; vMe(1) ; vMe(2) ; vMe(3); aMe(1) ;aMe(2);aMe(3);vV(1);vV(2);vV(3);aV(1);aV(2);aV(3);vE(1);vE(2);vE(3);aE(1);aE(2);aE(3);vM(1);vM(2);vM(3);aM(1);aM(2);aM(3);vJ(1);vJ(2);vJ(3);aJ(1);aJ(2);aJ(3)];