Extension to non-autonomous Hamiltonian systems

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 H) to account for time-dependence of the Hamiltonian:

 <* 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 include a line to increment the time value at each intermediary stage:

         t = t + c(i)*h

for i = 1(1)s. Finally a modified call to rksub.f is required

         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.