c
c     Numerical Analysis:
c     The Mathematics of Scientific Computing
c     D.R. Kincaid & E.W. Cheney
c     Brooks/Cole Publ., 1990
c
c     Section 8.3
c
c     Runge-Kutta-Fehlberg method 
c     for solving an initial value problem   
c
c
c     file: rkfelbergg.f
c
      external f
      dimension a(6),e(6),c(6),d(6,6),p(6)
      data M/500/, t/1.0/, x/2.0/, h/0.1/, delta/5.0e-5/, b/3.0/
c
      data (a(i),i=1,6)/
     + .118518518518518518518518519, 0.,  .518986354775828460038986355,
     + .50613149034201665780613149, -.18, .0363636363636363636363636364/
      data (c(i),i=1,6)/
     + 0., .25, .375, .923076923076923076923076923, 1., .5/
      data (e(i),i=1,6)/
     + .00277777777777777777777777778,0.,-.0299415204678362573099415205,
     + -.0291998936735778841041998937,.02,.0363636363636363636363636364/
      data ((d(i,j),j=1,6),i=1,6)/
     + 0., 0., 0., 0., 0., 0.,
     + .25, 0., 0., 0., 0., 0., 
     + .09375, 0., .28125, 0., 0., 0.,
     + .879380974055530268548020027,0.,-3.27719617660446062812926718,
     + 3.32089212562585343650432408, 0., 0.,
     + 2.03240740740740740740740741,0.,-8.,7.17348927875243664717348928,
     + -.205896686159844054580896686, .45297270955165692007797271,
     + -.296296296296296296296296296,0.,2.,
     + -1.38167641325536062378167641, 0., -.275/
c
      print *
      print *,' Runge-Kutta-Fehlberg method'
      print *,' Section 8.3, Kincaid-Cheney'
      print *
      print 7,'k','t','x','e','h' 
      iflag = 1 
c
      do 6 k=1,M
         xd = b-t
         if (abs(xd) .le. h) then
            iflag = 0
            h = xd
         endif
         y = x
         s = t
c
         do 3 i=1,6
            sum = 0.0
            do 2 j=1,i-1
               sum = sum + d(i,j)*p(j)
 2          continue
            p(i) = h*f(t+(c(i)*h),x+sum)
 3       continue
c
         sum = 0.0
         do 4 i=1,6
            sum = sum + a(i)*p(i)
 4       continue
         x = x + sum
c 
         est = 0.0
         do 5 i=1,6
            est = est + e(i)*p(i)
 5       continue
c
         t = t + h
         print 8,k,t,y,est,h 
c
         if (iflag .eq. 0) then
            stop
         endif
c
         if (abs(est) .ge. delta) then
            print *,' half h'
            h = h/2.0
            t = s
            x = y
            k = k-1
         else 
            if (abs(est) .lt. delta/128.0) then 
               print *,' double h'
               h = 2*h
            endif   
         endif
 6    continue
c
 7    format(a6,a9,2a15,a15)
 8    format(1x,i5,2x,4(e13.6,2x))
      stop
      end 
c  
      function f(t,x) 
      f = (t*x - x**2.0)/t**2.0
      return
      end