subroutine rkg(derive, x0, y, n, ndim, h, nx, f)
implicit real*8 (a-h, o-z)
external derive
dimension y(ndim, 0:*), f(ndim, 4)
data fp / 1. 70710 67811 86547 52440 d0 /
data fm / 0. 29289 32188 13452 47560 d0 /
data fn / 0. 16666 66666 66666 66667 d0 / ! 1/6
data fp2 / 3. 41421 35623 73095 04880 d0 / ! 2 + sqrt(2)
data fm2 / 0. 58578 64376 26904 95119 d0 /
data fpp / -4. 12132 03435 59642 57320 d0 /
data fmm / 0. 12132 03435 59642 57320 d0 /
hh = 0.5d0*h
fmh = fm*h
fph = fp*h
fnh = fn*h