***ADAMS-BASHFORTH********************************************************
*  PROGRAM FOR NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATION BY   *
* MULTI STEP ADAMS-BASHFORTH METHOD WITH RK-4 METHOD FOR STARTING VALUES *
**************************************************************************
*
      program ADAMBASH
      implicit doubleprecision(a-h,o-z)
      doubleprecision X(5),Y(5)
      external FXY
*
      print*, 'Enter initial and final values of independent variable
     $(Domain for which the solution is desired):'
      read(*,*) X(1), xu
      print*, 'Enter initial value of dependent variable:'
      read(*,*) Y(1)
      print*, 'Enter step size:'
      read(*,*) h
      n = (xu - X(1))/h
      write(6,*)'--------------------------------------------------'
      write(6,11) 'Independent Variable(x)','Dependent Variable(y)'
      write(6,*)'--------------------------------------------------'
      write(6,12) X(1),Y(1)
*
      do i = 1, 3
         CALL RK4(FXY,X(i),Y(i),h,slope)
         Y(i+1) = Y(i) + slope*h
         X(i+1) = X(i) + h
         write(6,12) X(i+1),Y(i+1)
      enddo 
*
      do i = 4, n
         CALL ADAMB4(FXY,X,Y,slope)
         Y(5) = Y(4) + slope*h
         X(5) = X(4) + h
         do k = 1, 4
            X(k) = X(k+1)
            Y(k) = Y(k+1)
         enddo
         write(6,12) X(4),Y(4)
      enddo
      write(6,*)'--------------------------------------------------'
*
11    format('',a22,a26)
12    format('',7x,f7.3,16x,f12.5)
*
      stop
      end
*
***RK-4th ORDER (CLASSIC METHOD)*************************************
*
      subroutine RK4(F,x,y,h,slope)
      implicit doubleprecision(a-h,o-z)
      doubleprecision k1,k2,k3,k4
*
      k1 = F(x,y)
      k2 = F(x+ h/2 , y+ h*k1/2)
      k3 = F(x+ h/2 , y+ h*k2/2)
      k4 = F(x+h , y+ h*k3)
      slope = (k1 + 2*k2 + 2*k3 + k4)/6d0
*
      return
      end
*
***4-th ORDER ADAMS-BASHFORTH******************************************
*
      subroutine ADAMB4(F,X,Y,slope)
      implicit doubleprecision(a-h,o-z)
      doubleprecision X(20),Y(20)
*
      slope = (55*F(X(4),Y(4)) - 59*F(X(3),Y(3)) + 37*F(X(2),Y(2)) - 9*F
     $(X(1),Y(1))) /24d0
*
      return
      end
*
***5-th ORDER ADAMS-BASHFORTH******************************************
*
      subroutine ADAMB5(F,X,Y,slope)
      implicit doubleprecision(a-h,o-z)
      doubleprecision X(20),Y(20)
*
      slope = (1901*F(X(5),Y(5)) - 2774*F(X(4),Y(4)) + 2616*F(X(3),Y(3))
     $- 1274*F(X(2),Y(2)) + 251*F(X(1),Y(1))) /720d0
*
      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-300d0)))
*
      return
      end