%Ode For Moon Orbits
function dYdt = moon_ode(t,Y)

r_io = [Y(1);Y(2)];
v_io = [Y(3);Y(4)];
r_eu = [Y(5);Y(6)];
v_eu = [Y(7);Y(8)];
r_gy = [Y(9);Y(10)];
v_gy = [Y(11);Y(12)];
r_j=[Y(13);Y(14)];
v_j=[Y(15);Y(16)];

% must declare them 'global' here and there
global Mj mio meuropa mganymede G 
%Forces "on what by what"
Fio_eu = (G*mio*meuropa)/(norm(r_eu-r_io)^3)*(r_eu-r_io);
Feu_io=-Fio_eu;
Fio_gy= (G*mio*mganymede)/(norm(r_gy-r_io)^3)*(r_gy-r_io);
Fgy_io=-Fio_gy;
Fio_j= (G*mio*Mj)/(norm(r_j-r_io)^3)*(r_j-r_io);
Fj_io=-Fio_j;
Feu_gy= (G*meuropa*mganymede)/(norm(r_gy-r_eu)^3)*(r_gy-r_eu);
Fgy_eu=-Feu_gy;
Feu_j= (G*meuropa*Mj)/(norm(r_j-r_eu)^3)*(r_j-r_eu);
Fj_eu=-Feu_j;
Fgy_j= (G*mganymede*Mj)/(norm(r_j-r_gy)^3)*(r_j-r_gy);
Fj_gy=-Fgy_j;

%accelerations
a_io=(Fio_eu+Fio_gy+Fio_j)/mio;
a_eu=(Feu_io+Feu_gy+Feu_j)/meuropa;
a_gy=(Fgy_io+Fgy_eu+Fgy_j)/mganymede;
a_j=(Fj_io+Fj_eu+Fj_gy)/Mj;


dYdt = [v_io(1),v_io(2),a_io(1),a_io(2),v_eu(1),v_eu(2),a_eu(1),a_eu(2) ...
    v_gy(1),v_gy(2),a_gy(1),a_gy(2),v_j(1),v_j(2),a_j(1),a_j(2)]';