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     Solving an initial value problem using Runge-Kutta method of order 4
c
c
c     file: rk4.f
c
      data M/256/, t/1.0/, x/2.0/, h/7.8125e-3/
c
      print *
      print *,' Runge-Kutta method of order 4'
      print *,' Section 8.3, Kincaid-Cheney'
      print *
      print 3,'k','t','x','e'
c
      e = abs(u(t) - x)
      print 4,0,t,x,e 
      h2 = 0.5*h
c
      do 2 k = 1,M
         f1 = h*f(t,x)       
         f2 = h*f(t + h2,x + 0.5*f1)   
         f3 = h*f(t + h2,x + 0.5*f2)   
         f4 = h*f(t + h,x + f3)
         x = x + (f1 + 2.0*(f2 + f3) + f4)/6.0
         t = t + h 
         e = abs(u(t) - x)
         print 4,k,t,x,e 
 2    continue
c
 3    format(a6,a9,2a16)
 4    format(1x,i5,2x,3(e14.7,2x))    
      stop
      end 
c
      function f(t,x) 
      f = (t*x -x**2.0)/t**2.0
      return
      end 
c  
      function u(t)
      u = t/(0.5 + alog(t))
      return
      end