Note: Orbits needs more accuracy near Moon

[rslippert]

Note: Current solve_ivp needs more accuracy near the moon,
where the derivatives are changing very fast
you will need to break up time evaluation into discrete t_points with t_eval=t_points

t_points = np.linspace(t_0, t_f, 4000)
sol = solve_ivp(EMBR0, [t_0, t_f], Y_0, atol=1e-13, rtol=1e-13,\
            method='LSODA', t_eval=t_points)