I can't figure out how to use this code for fixed step size.  imprk2/imprk4 seems to be not working at all for this case:
$ ./bug
Segmentation fault

Other steppers works (e.g. rk2) just fine.  See bug-ode2.c.

Second, Newton-type iterations could be inefficient in general.  See: Ernst Hairer, Christian Lubich, Gerhard Wanner, Geometric Numerical Integration.  Structure-Preserving Algorithms for Ordinary Differential Equations. Springer Series in Comput. Mathematics, Vol. 31, Springer-Verlag 2002. pp. 330-332.


(file #21314)

Quote from http://lists.gnu.org/archive/html/help-gsl/2009-12/msg00047.html

--->8---
It seems, for example, the 2 stages Gauss points are only iterate 3 steps, without checking whether it is converged or not.

The reason I care about this is, Implicit Gauss RK4 is symplectic, which for Hamiltonian Dynamics simulation is a very important property, it has better conservation of the first integral.
--->8---

Attached example code (bug.c) for Henon-Heiles problem shows the linear error grow in the Hamiltonian.  Compare with correct behavior for patched solvers (rk2imp.png/rk4imp.png).

Please, review attached patch, targeted to solve the problem.