************************************************************************
*                                                                      *
*     Subroutine GBS - Beta Version                                    *
*                                                                      *
*     GBS is Gragg-Bulirsch-Stoer integrator                           *
*     for solving first-order differential equations                   *
*     with variable step                                               *
*                                                                      *
*     Input parameters:                                                *
*     X     array of  integrated variables at TS                       *
*     TS    initial time                                               *
*     TF    final time                                                 *
*     HX    initial step size (autoselection if HX=0)                  *
*     NOR   order of integrator (NOR=4,6,8,...)                        *
*     NE    number of equations (integrated variables)                 *
*     ETOL  local error (constant HX if ETOL<=0)                       *
*     FUN   name of subroutine of right-hand members                   *
*           (EXTERNAL FUN required where GAUSS_32 called)              *
*                                                                      *
*     Output parameters:                                               *
*     X     array of integrated variables at TF                        *
*     HX    penultimate step size                                      *
*     NS    number of steps                                            *
*     NF    number of calls to subroutine FUN                          *
*                                                                      *
************************************************************************
      SUBROUTINE GBS(X,TS,TF,HX,NOR,NE,ETOL,NS,NF,FUN)
      ! Gragg-Bulirsch-Stoer Method
      ! Written by Avdyushev V. (2012)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(HMN=1D-8,SGM=3.16D0)
      DIMENSION NST(NOR/2),Y(NOR/2,NOR/2,NE)
      REAL*8, DIMENSION(NE) :: X,X0,X1,X2,F,F0
      LOGICAL LS,VS

      EXTERNAL FUN

      NF=0; NS=0

      IF(TS==TF) RETURN

      NC=NOR/2 ! Number of Table Columns

      IF(NC<=1) STOP 'Inadmissible Order!'

      ! Sequence of Spacings
      
      DO I=1,NC; NST(I)=2*I; END DO

      VS=ETOL>0D0 ! Variable Step Size
    
      ! Step for Zero Order
      
      IF(HX==0D0) THEN

      IF(.NOT.VS) STOP 'Zero Constant Step!'

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

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

      END IF

      T=TS; HX=DABS(HX)*(TF-TS)/DABS(TF-TS)
      DEG=1D0/DFLOAT(2*NC-1)
      RMX=SGM**DEG; RMN=1.D0/RMX; R=1D0
      LS=.FALSE.

      ! Integration

      DO

      CALL FUN(T,X,F0); NF=NF+1

    1 DO I=1,NC
      
      H=HX/DFLOAT(NST(I)); X0=X

      ! Euler Step

      X1=X0+H*F0; T1=T+H; J=0

      ! Midpoint Steps

      DO
      
      CALL FUN(T1,X1,F); NF=NF+1
      X2=X0+2D0*H*F; J=J+1

      IF(J==NST(I)) EXIT
      
      T1=T1+H; X0=X1; X1=X2

      END DO

      ! Smoothing Step

      Y(I,1,:)=(X0+2D0*X1+X2)/4D0

      END DO

      ! Aitken-Neville Algorithm
      
      DO J=1,NC-1
      DO I=J+1,NC
      DNM=(DFLOAT(NST(I))/DFLOAT(NST(I-J)))**2-1D0
      Y(I,J+1,:)=Y(I,J,:)+(Y(I,J,:)-Y(I-1,J,:))/DNM
      END DO
      END DO

      ! Step Control
      
      IF(VS) THEN
      
      ! Error Estimation
      
      ECAL=0D0

      DO I=1,NE
      ECAL=ECAL+(Y(NC,NC,I)-Y(NC,NC-1,I))**2
      END DO

      ECAL=DSQRT(ECAL)

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

      END IF

      ! Selection of First Step Size
      
      IF(NS==0.AND..NOT.LS) THEN
      IF(R>=RMX.OR.R<RMN) THEN
      HX=HX*R; GOTO 1
      END IF
      IF((TS-TF)*(TS-TF+HX)<=0D0) THEN
      HX_=HX; HX=TF-TS; LS=.TRUE.; GOTO 1
      END IF
      END IF

      ! Numerical Solution

      T=T+HX; X=Y(NC,NC,:); HX=HX*R; NS=NS+1

      IF(LS) EXIT

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

      HX=HX_

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