c=======================================================================
c     This is a subroutine that defined the equations of motion for
c     the driven harmonic oscillator system:
c
c     H = p**2/(2*m) + m*w0**2*x**2/2 + V*x*cos(w*t+theta)
c
c     The equations are set up to calculate the following quantities:
c
c     y(1) = S(x0,p0,t) (S is the action along the classical path)
c     y(2) = x(x0,p0,t)
c     y(3) = p(x0,p0,t)
c     y(4) = x(x0,p0+dp,t)
c     y(5) = p(x0,p0+dp,t)
c     y(6) = x(x0+dx,p0,t)
c     y(7) = p(x0+dx,p0,t)
c=======================================================================
      subroutine fcn(neq, t, y, yprime)

         implicit none

c-----------------------------------------------------------------------
c     Common block for communicating with orbitnt, jac
c-----------------------------------------------------------------------
         include 'fcn.inc'

c-----------------------------------------------------------------------
c     neq: number of odes (7)
c     t: time variable
c     y: solution of odes
c     yprime: derivative of solution
c-----------------------------------------------------------------------
         integer neq
         real*8 t
         real*8 y, yprime
         dimension y(1), yprime(1)

c-----------------------------------------------------------------------
c     Define derivatives
c-----------------------------------------------------------------------
         yprime(1) = y(3)**2/(2.0d0*m) - m*w0**2*y(2)**2 -
     &        V*y(2)*cos(w*t+theta)
         yprime(2) = y(3)/m
         yprime(3) = -m*w0**2*y(2) - V*cos(w*t+theta)
         yprime(4) = y(5)/m
         yprime(5) = -m*w0**2*y(4) - V*cos(w*t+theta)
         yprime(6) = y(7)/m
         yprime(7) = -m*w0**2*y(6) - V*cos(w*t+theta)
         
         return
      end


c=======================================================================
c     Optional Jacobian for lsoda
c
c     Note: in 'orbitnt.f' lsoda is set up for a banded Jacobian with
c     ml=1 and mu=2
c=======================================================================
      subroutine jac(neq,t,y,ml,mu,pd,nrowpd)
         
         implicit none

c-----------------------------------------------------------------------
c     Common communication with orbitnt and fcn
c-----------------------------------------------------------------------
         include 'fcn.inc'

c-----------------------------------------------------------------------
c     neq: number of odes (7)
c     ml: number of lower diagonals in Jacobian (1)
c     mu: number of upper diagonals in Jacobian (2)
c     nrowpd: number of rows in array pd (nrowpd=ml+mu+1=4)
c-----------------------------------------------------------------------
         integer neq
         integer ml, mu, nrowpd

c-----------------------------------------------------------------------
c     t: time variable
c     y: solution of odes
c     pd: array of partial derivatives (Jacobian)
c-----------------------------------------------------------------------
         real*8 t
         real*8 y, pd
         dimension y(1), pd(nrowpd,1)

c-----------------------------------------------------------------------
c     Define partial derivatives
c
c     note J(i,j) = pd(i-j+mu+1,j) (banded matrix form)
c-----------------------------------------------------------------------
         pd(2,2) = -2.0d0*m*w0**2*y(2) - V*cos(w*t+theta)
         pd(1,3) = y(3)/m
         pd(2,3) = 1.0d0/m
         pd(4,2) = -m*w0**2
         pd(2,5) = 1.0d0/m
         pd(4,4) = -m*w0**2
         pd(2,7) = 1.0d0/m
         pd(4,6) = -m*w0**2

         return
      end