      program RKInt
*---------------------------------------------------------------
*
*     Runge-Kutta-Integration 4. Ordnung
*
*---------------------------------------------------------------
      implicit none
      integer    Step,RKStep
      integer    DIM
      parameter  (DIM=4)
*
      real*8     vars(DIM)        ! p1,p2,x1,x2
c      real*8     dvars(DIM)       ! Ableitungen werden in DERIVS def.
      real*8     Time,dt,TimeStop
*
      character  Version*50
*---------------------------------------------------------------
      Version= '# RKInt 01 Juergen Schnack  02. 07. 2009'
*---------------------------------------------------------------
*
*     Anfangsbedingungen
*
      Time     =  0.d0
      vars(1)  =  0.d0           ! p1
      vars(2)  =  0.d0           ! p2
      vars(3)  = -2.d0           ! x1
      vars(4)  =  1.d0           ! x2

      dt       = 0.1d0
      TimeStop = 10.0d0
      Step     = 0
      RKStep   = 4

      write(6,'(a10,4a20)') 't','p1','p2','x1','x2'
      write(6,'(5f20.10)') Time,vars(1),vars(2),vars(3),vars(4)
*---------------------------------------------------------------
*
*     Zeitentwicklung
*
 200  continue
*
      Step = Step + 1
      Time = Time + dt

      call RKDUMB(vars,4,Time-dt,Time,RKStep) ! 4=DIM

      write(6,'(5f20.10)') Time,vars(1),vars(2),vars(3),vars(4)

      IF ((TimeStop-Time).gt.0) GOTO 200
*
 9999 continue
*
      stop
      end
***************************************************************
*                          ENDE VON MAIN
***************************************************************
*
*=== === === === === === === === === === === === === === === ==
*                          UNTERROUTINEN
*=== === === === === === === === === === === === === === === ==
*
*
*---------------------    DERIVS      ------------------------
*
      subroutine DERIVS(DIM,Time,vars,dvars)
*-------------------------------------------------------------
*
*     Derivs berechnet die ersten Ableitungen dvars
*     des Vektors vars nach t.
*
*---------------------------------------------------------------
      implicit none
      integer    DIM
      real*8     vars(DIM)        ! p1,p2,x1,x2
      real*8     dvars(DIM)       ! Ableitungen
      real*8     Time
*---------------------------------------------------------------
*     
*     Ableitungen (aus Mathematica)
*     
      dvars(1) = -1.d0 - vars(3) - 2.00485d0*(vars(3) - vars(4))
      dvars(2) = 2.00485d0*(vars(3) - vars(4)) - 4.d0*(-1.d0 + vars(4))
      dvars(3) = vars(1)
      dvars(4) = vars(2)
*     
*---------------------------------------------------------------
      return
      end
