I want to solve a relatively small system of stiff ODEs (< 10 first-order equations) using high precision floating point arithmetic (using MPFR or alike). What would be the easiest algorithm to implement/port? I'm not only interested in the "stepper" algorithm, but in the combination of stepper/error estimation/step size control.
I don't know if your problem would be tractable using a high order Taylor series method, but if so, and you're comfortable with implementing a solution to your problem in Python, then you could try a combination of mpmath, SymPy, and possibly pytaylor (which is not as well-established as the first two Python modules, but does implement 4th order RK in addition to Taylor series methods). Based on a quick read of a paper by Nedialkov in BIT, it seems as though a sufficiently high Taylor series (of order 20-30) could be used for highly accurate, efficient numerical integration of mildly to moderately stiff problems. (Nedialkov studies the differential-algebraic case; presumably, these results must also apply to ODEs.)
Whatever solution method you choose, I encourage you open-source your work, since it sounds like it'd be a useful contribution to the ODE software landscape.
This may be a bit outdated, but Hairer and Wanner's book recommends their own
radau5 code. In my own experience, this code is both extremely robust and efficient. It is also a rather straight-forward Runge-Kutta scheme, and thus not too difficult to implement.
There's a Matlab implementation of
radau5 by Christian Engstler somewhere out there.
Obviously, using multiple precision data types like those in MPFR only makes sense if numerical roundoff is at least around the same order of magnitude as the discretization error. You can only achieve this with high order integrators, so anything that has less than, say, 4th order will not be a useful candidate, and even higher order is better.
I'm not an expert in ODEs to recommend an integration package, but the considerations above may narrow your choice.
There are existing implementations, though I don't know of any that are publicly available. You could ask Glaser and Rohklin for their implementation of the methods described in this paper.
Is DLSODE an option?
Is the system linear? If so, DGPADM is probably the best.
They are relatively self-contained, so not too torturous to re-code if that's what you need to do.