Assume that the previous problem is extended to non-autonomous Hamiltonians. It is straightforward to modify the above subroutine for this case. In the template file rksub.mf, modify the first line to include time-dependence:
subroutine evalf(f,q,p,t)
The definition of f needs to be modified (the vector of first derivatives
of
<* FortranAssign[ f,
Flatten[{Outer[D,{H},pvars],- Outer[D,{H},qvars],D[H,t]}]
] *>
The dimension of f must also be increased:
f(<* 2 d + 1 *>)
The main routine should now includes a line to increment the time value at each
intermediary stage:
t = t + c(i)*h
for
call rksub(f,q,p,t)
To take the template file approach to an extreme, we could even select an implicit
or explicit Runge-Kutta method, depending upon some stiffness criteria
implemented in Mathematica.