***RUNGE-KUTTA-FEHLBERG**************************************************
*    PROGRAM FOR NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATION   *
*       RUNGE-KUTTA-FEHLBERG(5-th Order) METHOD WITH ERROR ESTIMATE     *
*************************************************************************
*
      program FEHLBERG
      implicit doubleprecision(a-h,o-z)
*
      external FXY
*
      print*, 'Enter initial and final values of independent variable
     $(Domain for which the'
      print*, 'solution is desired):'
      read(*,*) x, xu
      print*, 'Enter initial value of dependent variable:'
      read(*,*) y
      print*, 'Enter step size and print interval:'
      read(*,*) h, d
*
10    if (dabs(d) < dabs(h)) then
         write(6,*)
         print*, '*Print Interval Cannot be Less Than Step Size*'
         write(6,*)
         print*, 'Enter step size and print interval:'
         read(*,*) h, d
         goto 10
      endif
*
      m = (xu - x)/d
      n = d/h
      write(6,*)'-------------------------------------------------------
     $-----------------'
      write(6,101) 'Independent Variable(x)','Dependent Variable(y)','Er
     $ror Estimate'
      write(6,*)'-------------------------------------------------------
     $-----------------'
*
      write(6,102) x,y,0d0
      do i = 1, m
      do j = 1, n
         CALL FEHL(FXY,x,y,h,slope,error)
         y = y + slope*h
         x = x + h
      enddo
      write(6,102) x,y,error
      enddo
*
      write(6,*)'------------------------------------------------------'
101   format('',a23,a26,a21)
102   format('',7x,f7.3,16x,f12.4,16x,f8.5)
*
      stop
      end
*
**RUNGE-KUTTA-FEHLBERG(5-th ORDER)***********************************
*
      subroutine FEHL(F,x,y,h,slope,error)
      implicit doubleprecision(a-h,o-z)
      doubleprecision k1,k2,k3,k4,k5,k6
*
       k1 = F(x,y)
       k2 = F(x+ h/4d0 , y+ h*k1/4d0)
       k3 = F(x+ 3*h/8d0 , y+ h*(3*k1 + 9*k2)/32d0)
       k4 = F(x+ 12*h/13d0 , y+ h*(1932*k1 - 7200*k2 + 7296*k3)/2197d0)
       k5 = F(x+ h , y+ h*(439*k1/216d0 - 8*k2 + 3680*k3/513d0 - 
     $845*k4/4104d0))
       k6 = F(x+ h/2d0 , y+ h*(-8*k1/27d0 + 2*k2 - 3544*k3/2565d0 + 
     $1859*k4/4104d0 - 11*k5/40d0))
       slope = 16*k1/135d0 + 6656*k3/12825d0 + 28561*k4/56430d0 - 
     $9*k5/50d0 + 2*k6/55d0
       error = (k1/360d0 - 128*k3/4275d0 - 2197*k4/75240d0 + k5/50d0 + 
     $k6/27.5d0)* h
*
      return
      end
*
***THE DIFFERENTIAL EQUATION TO BE SOLVED; WRITTEN IN THE FORM:-
***                 dy/dx = F(x,y) = FXY
*
      function FXY(x,y)
      implicit doubleprecision(a-h,o-z)
*
      FXY = -1d0/(0.004d0* (1005.4d0 + 0.014999d0*(y-260d0) + 
     $0.00037498d0*(y-260d0)*(y-280d0) + 
     $0.00000d0*(y-260d0)*(y-280d0)*(y-300)))
*
      return
      end