************************************************************************
*                                                                      *
*     Subroutine AMB - Beta Version                                    *
*                                                                      *
*     AMB is Adams-Moulton-Bashforth integrator                        *
*     for solving first-order differential equations                   *
*     with variable step                                               *
*                                                                      *
*     Input parameters:                                                *
*     TS    initial time                                               *
*     TF    final time                                                 *
*     X     array of  integrated variables at TS                       *
*     NOR   order of integrator (unrestricted)                         *
*     NE    number of equations (integrated variables)                 *
*     ETOL  local error                                                *
*     NI    number of iterations for corrector                         *
*     FUN   name of subroutine of right-hand members                   *
*           (EXTERNAL FUN required where AMB called)                   *
*                                                                      *
*     Output parameters:                                               *
*     X     array of integrated variables at TF                        *
*     NST   number of steps                                            *
*     NF    number of calls to subroutine FUN                          *
*                                                                      *
************************************************************************
      SUBROUTINE AMB(X,TS,TF,NOR,NE,ETOL,NI,NST,NF,FUN)
      ! Adams-Moulton-Bashforth Method
      ! Written by Avdyushev V. (2012)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(SGM=3.16D0,HMN=1D-8)
      DIMENSION B(0:NOR-1),G(NOR),D(0:NOR-1,NE)
      DIMENSION X(NE),XP(NE),XC(NE),F(0:NOR-1,NE),T(0:NOR-1)
      LOGICAL LS
 
      EXTERNAL FUN

      NF=0; NST=0

      IF(TS==TF) RETURN

      NS=NOR-1
      IF(NS<1) STOP 'Inadmissible Order!'
      IF(NI<1) NI=1
      ETOL=DABS(ETOL)

      ! Step for Zero Order

      CALL FUN(TS,X,F(0,:)); NF=1

      ECAL=0D0
      DO I=1,NE
      ECAL=ECAL+F(0,I)**2
      END DO
      IF(ECAL/=0D0) THEN
      H=ETOL/DSQRT(ECAL)
      ELSE
      H=HMN
      END IF

      H=H*(TF-TS)/DABS(TF-TS)

      ! Initialization

      R=1D0; N=1; NPO=2; NMO=0
      DEG=0.5D0; RMX=SGM**DEG; RMN=1.D0/RMX
      T(0)=TS; LS=.FALSE. 
    1 H=H*R 
    2 T(1)=TS+H; XC=X

      DO ! Integration

      ! Gammas

      TMP=1D0
      DO J=1,NPO
      TMP=TMP*H/DFLOAT(J)
      G(J)=TMP
      END DO
      B(0)=G(1)

      DO I=1,N
      DO J=1,NPO-I
      G(J)=(T(N)-T(N-I))*G(J)-DFLOAT(J)*G(J+1)
      END DO
      B(I)=G(1)
      END DO

      ! Bashforth (Predictor)

      DO I=0,NMO
      D(I,:)=F(NMO-I,:)
      DO J=0,I-1
      TMP=T(NMO-I)-T(NMO-J)
      D(I,:)=(D(I,:)-D(J,:))/TMP
      END DO
      END DO

      DO I=0,NMO
      XC=XC+B(I)*D(I,:)
      END DO

      XP=XC ! Predictor

      DO K=1,NI ! Iterations

      ! Force at Point N

      CALL FUN(T(N),XC,F(N,:)); NF=NF+1

      ! Moulton (Corrector)

      D(N,:)=F(N,:)
      DO J=0,NMO
      TMP=T(N)-T(NMO-J)
      D(N,:)=(D(N,:)-D(J,:))/TMP
      END DO

      XC=XP+B(N)*D(N,:) ! Corrector

      END DO

      IF(LS) EXIT

      ! Step Control

      ECAL=0D0
      DO I=1,NE
      ECAL=ECAL+D(N,I)**2
      END DO
      ECAL=DABS(B(N))*DSQRT(ECAL)

      IF(ECAL/=0D0) THEN
      R=(ETOL/ECAL)**DEG
      IF(R>RMX) R=RMX
      ELSE
      R=RMX
      END IF

      IF(NST==0) THEN
      IF(R>=RMX.OR.R<RMN) GOTO 1 
      IF((TS-TF)*(TS-TF+H)<=0D0) THEN
      H=TF-TS; LS=.TRUE. ! Last Step
      GOTO 2
      END IF
      END IF

      IF(R>RMX) R=RMX; H=H*R; NST=NST+1

      IF((T(N)-TF)*(T(N)-TF+H)<=0D0) THEN
      H=TF-T(N); LS=.TRUE. ! Last Step
      END IF

      IF(N<NS) THEN ! Assembling of Integration Scheme 
      NMO=N; N=N+1; NPO=N+1; T(N)=T(NMO)+H
      DEG=1.D0/DFLOAT(NPO); RMX=SGM**DEG; RMN=1.D0/RMX
      ELSE ! Step-by-Step Advancement
      F(0:NMO,:)=F(1:N,:)
      T(0:NMO)=T(1:N); T(N)=T(N)+H
      END IF

      END DO

      X=XC; NST=NST+1 ! Final Result

      END
***********************************************************************
      SUBROUTINE FUN(T,X,F)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION X(6),F(6)
      R=DSQRT(X(1)**2+X(2)**2+X(3)**2)
      F(1:3)= X(4:6)
      F(4:6)=-X(1:3)/R**3
      END SUBROUTINE FUN
