!********************************************************************
!*    Solve Y' = F(X,Y) with Initial Condition Y(X0)=Y0 using the   *
!*    Adams-Bashforth Method                                        *
!* ---------------------------------------------------------------- *
!* REFERENCE: "Méthode de calcul numérique- Tome 2 - Programmes en  *
!*             Basic et en Pascal By Claude Nowakowski, Edition du  *
!*             P.S.I., 1984" [BIBLI 04].                            * 
!* ---------------------------------------------------------------- *
!* SAMPLE RUN:                                                      *
!* (Solve Y' = -Y + X/((1+X)*(1+X))  with Y(0)=1 for X0=0 up to     *
!*  X1=1.0, exact solution is Y = 1 / (1+X).                        *
!*                                                                  *
!*      X           Y        Y exact      Error                     *
!* -----------------------------------------------                  *
!*   0.000000    1.000000    1.000000   0                           *
!*   0.050000    0.952381    0.952381   0.000000                    *
!*   0.100000    0.909091    0.909091   0.000000                    *
!*   0.150000    0.869525    0.869565   0.000040                    *
!*   0.200000    0.833265    0.833333   0.000068                    *    
!*   0.250000    0.799910    0.800000   0.000090                    *
!*   0.300000    0.769125    0.769231   0.000106                    *
!*   0.350000    0.740623    0.740741   0.000117                    *
!*   0.400000    0.714160    0.714286   0.000125                    *
!*   0.450000    0.689525    0.689655   0.000131                    *
!*   0.500000    0.666533    0.666667   0.000134                    *
!*   0.550000    0.645026    0.645161   0.000135                    *
!*   0.600000    0.624865    0.625000   0.000135                    *
!*   0.650000    0.605926    0.606061   0.000134                    *
!*   0.700000    0.588103    0.588235   0.000133                    *
!*   0.750000    0.571298    0.571429   0.000131                    *
!*   0.800000    0.555428    0.555556   0.000128                    *
!*   0.850000    0.540416    0.540541   0.000125                    *
!*   0.900000    0.526194    0.526316   0.000121                    *
!*   0.950000    0.512703    0.512821   0.000118                    *
!*   1.000000    0.499886    0.500000   0.000114                    *
!* -----------------------------------------------                  *
!*                                                                  * 
!*                                F90 Release By J-P Moreau, Paris. *
!*                                        (www.jpmoreau.fr)         *
!********************************************************************
Program AdamsBashforth

implicit real*8 a-h,o-z

  dimension B(0:3),X(0:3),Y(0:3)

  H = 0.05d0  !integration step

  !Initial conditions
  X(0)=0.d0             !starting abscissa X0 
  X1=1.d0               !ending abscissa X1
  Y(0)=1.d0             !Initial value of Y at X0
  !write header
  print *,'       X           Y        Y exact      Error   '
  print *,'  -----------------------------------------------'
  !write initial line
  write(*,10) X(0), Y(0), Y(0),'    0'
  !use Runge-Kutta to start
  do K=0, 1
    call RK4(K,H,X,Y); YEX=FX(X(K+1)); ER=DABS(YEX-Y(K+1))
    write(*,11) X(K+1), Y(K+1), YEX, ER
  end do
  !main integration loop
  do while(X(2)<X1)
    do I=1, 3 
	  B(I)=F(X(3-I),Y(3-I))
    end do
    X(3)=X(2)+H
    Y(3)=Y(2)+H*(23.d0*B(1)-16.d0*B(2)+5.d0*B(3))/12.d0        
    YEX=FX(X(3)); ER=DABS(Y(3)-YEX)
    write(*,11) X(3), Y(3), YEX, ER
    do K=0, 2
      X(K)=X(K+1); Y(K)=Y(K+1)
    end do    
  end do
  print *,'  -----------------------------------------------'

10 format(F12.6,F12.6,F12.6,A5)
11 format(F12.6,F12.6,F12.6,F12.6)

END

    !User defined function Y'=F(X,Y)
    real*8 Function F(X,Y)
      real*8 X, Y
      F = -Y + X/((1.d0+X)*(1.d0+X))
    End

    !Exact solution Y=FX(X)
    real*8 Function FX(X)
      real*8 X
	  FX = 1.d0/(1.d0+X)
    End

    !Runge-Kutta method to calculate first two points only
    Subroutine RK4(K,H,X,Y)
	  real*8 H,X(0:3),Y(0:3)
	  real*8 C1,C2,C3,C4, F
      C1=F(X(K),Y(K))
      C2=F(X(K)+H/2.d0,Y(K)+H/2.d0*C1)
      C3=F(X(K)+H/2.d0,Y(K)+H/2.d0*C2)
      C4=F(X(K)+H,Y(K)+H*C3)
      X(K+1)=X(K)+H
      Y(K+1)=Y(K)+H*(C1+2*C2+2*C3+C4)/6.d0
	  Return
    End

!end of file adambash.f90