I would be very interested in why it would not work with Boost.Units since I might be interested in using such a library with my own datatypes instead of double or std::complex<double>
If you want to use Boost.Units the system function has to work on different types, for example: harmonic_oscillator( StateType , DerivType , TimeType ) The units have to obey [DerivType] = [StateType] / [TimeType]. This condition is not a problem for the explicit Runge-Kutta solvers. But if you want to use some implicit methods you are usually faced with some matrix inversion routines. In this case the state type has the restrictions to work with appropriate LA functions. I doubt that it will be possible to get this working with Boost.Units. Another problem might be that every entry in the state vector could have a different dimension, although this might be solved by using fusion vectors or something similar. If there is a strong interest in Boost.Units we could try to include them for some specific steppers in the next version. But in my opinion this is not that important. If you want to use your own datatypes you can try to specialize the container_traits. At the moment we are working on a small redesign of the library allowing for an stronger independence of the underlying types. One goal is for example to get the library working with thrust to enable CUDA support (http://code.google.com/p/thrust/). What kind of datatypes do you have in mind? Maybe we can try to get this library working with your types, too. Best regards, Karsten