implicit real ( a-h, o-z )
      real y(2 )
      external simple
      xi = 10
      xf = 0
      y(1) = exp(xi)
      y(2) = exp(xi)
      n = 2
      del = .1
      accurc = 1.e-8
      imax = 10
      write(*,*) kutta( xi, xf, y, n, del, accurc, imax, simple ),
     *              ' iterations'
      do 12500 i = 1, 2
         write(*,100) i, xf, y(i)-1
12500 continue
      stop
100   format( ' error in y(', i2, ') at x=', g12.5, ' is ', g12.5 )
      end
      subroutine simple( x, y, f )
      implicit real ( a-h, o-z )
      real y(2 ), f(2 )
      f(1) = y(2)
      f(2) = 2*y(2) - y(1)
      return
      end
      integer function kutta( xi, xf, y, n, del, accurc, imax, equa )
      implicit real (a-h,o-z)
      real   y(6 ), yi(6 ), yn(6 ),
     *              k1(6 ), k2(6 ), k3(6 ), k4(6 ), k5(6 ),
     *              f(6 ), e(6 ), f1(6 )
      real   test, xi, xf, del, accurc, xn, h
      real   amax1, abs
      logical quit
      if ( n .gt.  6  ) then
         write(*,*) 'KUTTA: too many equations: ', n, '>', 6
         kutta = -1
         return
      end if
      kutta = 0
      nhalves = 0
      xn = xi
      h = del
      quit = .false.
      do 10000 i000 = 1, ( n ), 1
      yn(i000) = y(i000)
10000 continue
      if (  ( xf .gt.  xi .and.  h .lt.  0 )  .or.   ( xf .lt.  xi .and.
     *  h .gt.  0 )  ) h = -h
17500 continue
         if ( ( ( (xn+h.ge. xf) .and.  (xf.ge. xi) )  .or.  ( (xn+h.le. 
     *xf) .and.  (xf.le. xi) )  )  ) then
            del = h
            h = xf - xn
            quit = .true.
         end if
         call equa(xn, yn, f1)
17501 continue
      do 10001 i000 = 1, ( n ), 1
            k1(i000) = h*f1(i000)/3.
            yi(i000) = yn(i000) + k1(i000)
10001 continue
            call equa(xn + h/3., yi, f)
      do 10002 i000 = 1, ( n ), 1
            k2(i000) = h*f(i000)/3.
            yi(i000) = yn(i000) + k1(i000)/2. + k2(i000)/2.
10002 continue
            call equa(xn + h/3., yi, f)
      do 10003 i000 = 1, ( n ), 1
            k3(i000) = h*f(i000)/3.
            yi(i000) = yn(i000) + 3.*k1(i000)/8. + 9.*k3(i000)/8.
10003 continue
            call equa(xn + h/2., yi, f)
      do 10004 i000 = 1, ( n ), 1
            k4(i000) = h*f(i000)/3.
            yi(i000) = yn(i000) + 3.*k1(i000)/2. - 9.*k3(i000)/2. + 6.*k
     *4(i000)
10004 continue
            call equa(xn + h, yi, f)
            test = 0.0
      do 10005 i000 = 1, ( n ), 1
            k5(i000) = h*f(i000)/3.
            e(i000) = (k1(i000) - 9.*k3(i000)/2. + 4.*k4(i000) - k5(i000
     *)/2.)/5.
            test = amax1(test, abs(e(i000)))
10005 continue
      if (.not.( test .ge.  accurc )) goto 15001
            nhalves = nhalves + 1
            if( nhalves .ge.  imax ) then
               del = h
               kutta = -1
               write(*,*) 'KUTTA: step made too small:'
               write(*,*) '       del = ', del, '    at x = ', xn
               return
            end if
            h = h/2.
            quit = .false.
      goto 17501
15001 continue
      do 10006 i000 = 1, ( n ), 1
         yn(i000) = yn(i000) + (k1(i000) + 4.*k4(i000) + k5(i000))/2.
10006 continue
         xn = xn + h
         kutta = kutta + 1
         if ( test .lt.  accurc/32. ) h = 2.*h
      if (.not.( quit )) goto 17500
15000 continue
      do 10007 i000 = 1, ( n ), 1
      y(i000) = yn(i000)
10007 continue
      return
      end