Consider custom gradients for N-body problems

[pnkraemer]

Many differential equations, for instance the three-body problem or the Pleiades problem involve terms like

f(x, y) = x^1.5 / (x^2 + y^2)^0.5

which leads to differentiability issues around zero especially when considering high-order derivatives of the vector field. For this reason, e.g. https://pnkraemer.github.io/probdiffeq/benchmarks/pleiades/external/ defines its own Pleiades implementation which works with high-order Taylor series (the current implementation in DiffEqZoo does not work in this notebook).

One issue with simply copying the implementation in this notebook is that the numpy backend is not happy with dividing by zero.