      subroutine dopri5(n,t,tend,h,x,dgl,work,lwk,ierr)
C
C Integrator mit Dormand-Prince Koeffizienten 5ter Ordnung.
C Integriert das System DGL der Ordnung N von (T0, X) bis (TEND, X)
C mit Schrittweite H.
C
C     n.................... Systemdimension             (I)
C     t.................... aktuelle Zeit               (I/O)
C     tend................. Endzeit                     (I)
C     h.................... aktuelle Schrittweite       (I/O)
C     x.................... aktueller Zustand           (I/O)
C     dgl.................. Differentialgleichung       (I)
C     work................. Arbeitsfeld der Laenge lwk  (I/O)
C     lwk.................. Laenge des Arbeitsfelds     (I)
C     ierr................. Fehlerflag                  (O)
C                           =  0  success
C                           = -1  lwk zu klein
C                           = -2  h wird zu klein
C                           = -4  h0 zu klein (<= 0)
C
      implicit none
      external dgl
      integer n, lwk, ierr
      double precision t, tend, h, x(n), work(lwk)
C
C...Koeffizienten 5/4ter Ordnung nach Dormand und Prince
C
      double precision a21, a31, a32, a41, a42, a43, a51, a52, a53,
     &                 a54, a61, a62, a63, a64, a65, a71, a72, a73,
     &                 a74, a75, a76, b1hat, b2hat, b3hat, b4hat, 
     &                 b5hat, b6hat, b7hat, c2, c3,c 4, c5, c6, c7
      parameter (a21=1.d0/5.d0, a31=3.d0/40.d0 , a32=9.d0/40.d0, 
     &           a41=44.d0/45.d0, a42=-56.d0/15.d0, a43=32.d0/9.d0, 
     &           a51=19372.d0/6561.d0, a52=-25360.d0/2187.d0,
     &           a53=64448.d0/6561.d0, a54=-212.d0/729.d0, 
     &           a61=9017.d0/3168.d0, a62=-355.d0/33.d0, 
     &           a63=46732.d0/5247.d0, a64=49.d0/176.d0, 
     &           a65=-5103.d0/18656.d0, a71=35.d0/384.d0, a72=0.d0, 
     &           a73=500.d0/1113.d0, a74=125.d0/192.d0, 
     &           a75=-2187.d0/6784.d0, a76=11.d0/84.d0)
      parameter (b1hat=5179.d0/57600.d0, b2hat=0.d0,
     &           b3hat=7571.d0/16695.d0, b4hat=393.d0/640.d0,
     &           b5hat=-92097.d0/339200.d0, b6hat=187.d0/2100.d0,
     &           b7hat=1.d0/40.d0)
      parameter (c2=1.d0/5.d0, c3=3.d0/10.d0, c4=4.d0/5.d0, 
     &           c5=8.d0/9.d0, c6=1.d0, c7=1.d0)
C
C...Verschiedenes
      double precision tj, hact
      integer i
      logical last
C
C...Adresszeiger der sieben RK-Stufen, des Hilfszustands
C   und Ueberpruefung von lwk, h
C
      integer f1, f2, f3, f4, f5, f6, f7, xhatj
C
      f1 = 1
      f2 = f1 + n
      f3 = f2 + n
      f4 = f3 + n
      f5 = f4 + n
      f6 = f5 + n
      f7 = f6 + n
      xhatj = f7 + n
      if ((xhatj+n-1) .gt. lwk) then
         ierr = -1
         goto 999
      endif
      if (h.le.0.d0) then
         ierr = -4
         goto 999
      endif
C
C...Initialisierung
C
      ierr = 0
C
C...Ergebnisausgabe, waere bei richtigem Integrator nicht da!
C
      write(*,1111) t, (x(i), i = 1,n)
C
C...Entry point fuer den Integrationsloop
C
 100  continue
      if ((t+h) .ge. tend) then
         hact = tend-t
         last = .true.
      else
         hact = h
         last = .false.
      endif
C
C...Die sieben RK-Stufen
C
C -> f1
      do i=1,n
         work(xhatj+i-1) = x(i)
      enddo
      tj = t
      call dgl(n, tj, work(xhatj), work(f1))
C -> f2
      do i=1,n
         work(xhatj+i-1) = x(i)+hact*( a21*work(f1+i-1))
      enddo
      tj = t+c2*hact
      call dgl(n, tj, work(xhatj), work(f2))
C -> f3
      do i=1,n
         work(xhatj+i-1) = x(i)+hact*( a31*work(f1+i-1)
     &                                +a32*work(f2+i-1))
      enddo
      tj = t+c3*hact
      call dgl(n, tj, work(xhatj), work(f3))
C -> f4
      do i=1,n
         work(xhatj+i-1) = x(i)+hact*( a41*work(f1+i-1)
     &                                +a42*work(f2+i-1)
     &                                +a43*work(f3+i-1))
      enddo
      tj = t+c4*hact
      call dgl(n, tj, work(xhatj), work(f4))
C -> f5
      do i=1,n
         work(xhatj+i-1) = x(i)+hact*( a51*work(f1+i-1)
     &                                +a52*work(f2+i-1)
     &                                +a53*work(f3+i-1)
     &                                +a54*work(f4+i-1))
      enddo
      tj = t+c5*hact
      call dgl(n, tj, work(xhatj), work(f5))
C -> f6
      do i=1,n
         work(xhatj+i-1) = x(i)+hact*( a61*work(f1+i-1)
     &                                +a62*work(f2+i-1)
     &                                +a63*work(f3+i-1)
     &                                +a64*work(f4+i-1)
     &                                +a65*work(f5+i-1))
      enddo
      tj = t+c6*hact
      call dgl(n, tj, work(xhatj), work(f6))
C -> f7
      do i=1,n
         work(xhatj+i-1) = x(i)+hact*( a71*work(f1+i-1)
     &                                +a72*work(f2+i-1)
     &                                +a73*work(f3+i-1)
     &                                +a74*work(f4+i-1)
     &                                +a75*work(f5+i-1)
     &                                +a76*work(f6+i-1))
      enddo
      tj = t+c7*hact
      call dgl(n, tj, work(xhatj), work(f7))
C
C...Zusammenbau der neuen Loesung und Zeitfortschaltung
C
      do i=1,n
         x(i) = x(i)+hact*( b1*work(f1+i-1)
     &                     +b2*work(f2+i-1)
     &                     +b3*work(f3+i-1)
     &                     +b4*work(f4+i-1)
     &                     +b5*work(f5+i-1)
     &                     +b6*work(f6+i-1)
     &                     +b7*work(f7+i-1))
      enddo
      t = t+hact
C
C...Ergebnisausgabe, waere bei richtigem Integrator nicht da!
C
      write(*,1111) t, (x(i), i = 1,n)
C
C...Naechster Schritt
C
      if (.not.last) goto 100
C
C...Ende
C
 999  return
 1111 format (3E13.5)
      end
