function dYdt = kepler_ode_final(Y)
% assume 'Y' is a vector of form y = [Xh;Yh;VXh;VYh;Xe;Ye;Xm;Ym]

Rh = [Y(1);Y(2)];
Vh = [Y(3);Y(4)];
Re = [Y(5);Y(6)];
Rm = [Y(7);Y(8)];

global G Ms Me Mm Mh Feh Fsh

Fsh = -Ms * Mh * G / norm(Rh)^3 .* Rh;
Feh = -Me * Mh * G / norm(Rh-Re)^3 .* (Rh-Re);
Fmh = -Mm * Mh * G / norm(Rh-Rm)^3 .* (Rh-Rm);

Ah = (Fsh + Feh + Fmh) ./ Mh;

dYdt = [Vh(1);Vh(2);Ah(1);Ah(2);Re(1);Re(2);Rm(1);Rm(2)];