***HEUN'S METHOD********************************************************
*   PROGRAM FOR NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATION   *
*                BY HEUN'S PREDICTOR- CORRECTOR METHOD                 *
************************************************************************
*
      program HEUNS
      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 (abs(d) < abs(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
      print*, 'Maximum limit of iteration:'
      read(*,*) maxim
*
      m = (xu - x)/d
      n = d/h   
      write(6,*)'-------------------------------------------------'
      write(6,*) '    Indpndt variable(x)        Dpndt variable(y)'
      write(6,*)'-------------------------------------------------'
*
      write(6,101) x,y
      do i = 1, m
      do j = 1, n
         CALL HEUN(FXY,x,y,h,maxim,slope)
         y = y + slope*h
         x = x + h
      enddo
      write(6,101) x,y
      enddo
*
      write(6,*)'--------------------------------------------'
101   format('',8x,f7.2,14x,f11.4) 
*
      stop
      end
* 
***SUBROUTINE HEUN'S METHOD*******************************************
*
      subroutine HEUN(F,x,y,h,maxim,slope)
      implicit doubleprecision(a-h,o-z)
*
      slope1 = F(x,y)
      x1 = x + h
      y1 = y + slope1 *h
*
      do i = 1, maxim
         slope2 = F(x1,y1)
         slope  = (slope1 + slope2)/2d0
         y2 = y + slope *h
         y1 = y2
      enddo
*
      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