Has anyone here tried implementing the embedded 3(2) Runge-Kutta pair of Bogacki and Shampine in C or Fortran?

If not, how might I modify the NR routine rkqs() to be able to use this formula?

(For more info. on this, check the paper by Bogacki and Shampine).

Thanks a lot.

Jan M. (^_^)