ODE with externally calculated J^{-1}

[Eric Switzer]

I'm working on a problem that involves 10^6 coupled, stiff
differential equations.  The Jacobian inverse for the system is too
big to do within GSL, so I'd like to write a code that calculates
J^{-1} x using some tricks, instead of the LU method in bsimp.  The
options seem to be:
1. modify the GSL stepper and recompile GSL.  This makes the software
hard to distribute.
2. rewrite the ODE solver from scratch, (maybe more complicated than
NR, but less general than GSL?)
3. copy the GSL code to the radiative transfer directory, modify it,
and rename all the objects so that it doesn't clash with GSL? 
Something with pointers?

Is it true that you can't add additional methods without changing the
base GSL code?  Thanks,

Eric