这可能有点过时，但 Hairer 和 Wanner 的书推荐了他们自己的radau5代码。以我自己的经验，这段代码非常健壮和高效。这也是一个相当直接的龙格-库塔方案，因此实施起来并不难。

Christian Engstler 在某个地方有一个 Matlab 实现radau5。

显然，仅当数值舍入至少与离散化误差在同一数量级左右时，使用 MPFR 中的多精度数据类型才有意义。您只能通过高阶积分器来实现这一点，因此任何小于 4 阶的东西都不是有用的候选者，甚至更高阶更好。

我不是推荐集成包的 ODE 专家，但上述考虑可能会缩小您的选择范围。

有现有的实现，尽管我不知道任何公开可用的。您可以向 Glaser 和 Rohklin 询问他们对本文所述方法的实施情况。

我不知道使用高阶泰勒级数方法是否可以解决您的问题，但如果是这样，并且您对在 Python 中实现问题的解决方案感到满意，那么您可以尝试结合使用mpmath、SymPy和可能是pytaylor（它不像前两个 Python 模块那么完善，但除了泰勒级数方法之外，它确实实现了 4 阶 RK）。根据对Nedialkov 在 BIT 中的一篇论文的快速阅读，似乎足够高的泰勒级数（20-30 阶）可用于高度准确、有效的轻度到中度刚性问题的数值积分。（Nedialkov 研究微分代数情况；推测这些结果也必须适用于 ODE。）

无论您选择哪种解决方案方法，我都鼓励您开源您的工作，因为这听起来像是对 ODE 软件领域的有益贡献。