Hi there, I've been playing around with your very nice solver--thanks for putting those into the public domain!
I noticed that there's a sign error for the Jupiter force calculation. It probably doesn't matter for Jupiter mass planets, but when integrating planets closer to each other in mass, then that error will become important. The force function should be
Also, it would be beneficial to forward-integrate the motions of both planets in the Runge-Kutta solver symmetrically for cases where the planets are closer to each other and/or similar in mass. Here's a way to do that:
r1e=re+k11e*h*0.5
v1e=ve+k21e*h*0.5
r1j=rj+k11j*h*0.5
v1j=vj+k21j*h*0.5
k12e = dr_dt(t + 0.5*h,r1e,v1e,'earth',r1j,v1j)
k22e = dv_dt(t + 0.5*h,r1e,v1e,'earth',r1j,v1j)
k12j = dr_dt(t + 0.5*h,r1j,v1j,'jupiter',r1e,v1e)
k22j = dv_dt(t + 0.5*h,r1j,v1j,'jupiter',r1e,v1e)
r2e=re+k12e*h*0.5
v2e=ve+k22e*h*0.5
r2j=rj+k12j*h*0.5
v2j=vj+k22j*h*0.5
k13e = dr_dt(t + 0.5*h,r2e,v2e,'earth',r2j,v2j)
k23e = dv_dt(t + 0.5*h,r2e,v2e,'earth',r2j,v2j)
k13j = dr_dt(t + 0.5*h,r2j,v2j,'jupiter',r2e,v2e)
k23j = dv_dt(t + 0.5*h,r2j,v2j,'jupiter',r2e,v2e)
r3e=re+k13e*h
v3e=ve+k23e*h
r3j=rj+k13j*h
v3j=vj+k23j*h
k14e = dr_dt(t + h,r3e,v3e,'earth',r3j,v3j)
k24e = dv_dt(t + h,r3e,v3e,'earth',r3j,v3j)
k14j = dr_dt(t + h,r3j,v3j,'earth',r3e,v3e)
k24j = dv_dt(t + h,r3j,v3j,'earth',r3e,v3e)
y0e = re + h * (k11e + 2.*k12e + 2.*k13e + k14e) / 6.
y1e = ve + h * (k21e + 2.*k22e + 2.*k23e + k24e) / 6.
y0j = rj + h * (k11j + 2.*k12j + 2.*k13j + k14j) / 6.
y1j = vj + h * (k21j + 2.*k22j + 2.*k23j + k24j) / 6.
ze = np.zeros([2,2])
ze = [y0e, y1e]
zj = np.zeros([2,2])
zj = [y0j, y1j]
return ze,zj:
Finally, you can avoid the case distinction in the gravity term if you write the force terms as (e.g.):