Bulirsch-Stoer Step

SUBROUTINE bsstep(y, dydx, x, htry, eps, yscal, hdid, hnext, derivs)

  USE nrtype; USE nrutil, ONLY : arth, assert_eq, cumsum, iminloc, nrerror, &
       outerdiff, outerprod, upper_triangle
  USE nr, ONLY : mmid, pzextr
  IMPLICIT NONE
  REAL(sp), DIMENSION(:), INTENT(inout) :: y
  REAL(sp), DIMENSION(:), INTENT(in) :: dydx, yscal
  REAL(sp), INTENT(inout) :: x
  REAL(sp), INTENT(in) :: htry, eps
  REAL(sp), INTENT(out) :: hdid, hnext
  INTERFACE
     SUBROUTINE derivs(x, y, dydx)
       USE nrtype
       IMPLICIT NONE
       real9sp), INTENT(in) :: x
       REAL(sp), DIMENSION(:), INTENT(in) :: y
       REAL(sp), DIMENSION(:), INTENT(out) :: dydx
     END SUBROUTINE derivs
  END INTERFACE
  INTEGER(i4b), PARAMETER :: imax = 9, kmaxx=imax-1
  REAL(sp), PARAMETER :: safe1=0.25_sp, safe2=0.7_sp, redmax=1.0e-5_sp, &
       redmin=0.7_sp, tiny=1.0e-30_sp, scalmx=0.1_sp

  ! Bulirsch-Stoer step with monitoring of local truncation error to 
  ! ensure accuracy and adjust stepsize. Input are athe dependent variable
  ! vector y and its derivative dydx at the starting value of the independent
  ! variable x. Also input are the stepsize that was actually accomplished,
  ! and hnext is the estimated next stepsize. derivs is the user-supplied
  ! subroutine that computes the right-hand-side derivatives. y, dydx, and
  ! yscal must all have the same length. Be sure to set htry on successive 
  ! steps to the value of hnext returned from the previous step, as is the 
  ! case if the routine is called by odeint.
  ! Parameters: kmaxx is the maximum row number used in the extrapolation;
  ! imax is the next row number; safe1 and safe2 are safety factors;
  ! redmax is the maximum factor used when a stepsize is reduced, redmin
  ! the minimum; tiny prevents division by zero; 1/scalmx is the maximum
  ! factor by which a stepsize can be increased.

  INTEGER(i4b) :: k, km, ndum
  INTEGER(i4b, DIMENSION(imax) :: nseq = (/ 2, 4, 6, 8, 10, 12, 14, 16, 18 /)
  INTEGER(i4b), SAVE :: kopt, kmax
  REAL(sp), DIMENSION(kmaxx, kmaxx), SAVE :: alf
  REAL(sp), DIMENSION(kmaxx) :: err
  REAL(sp), DIMENSION(imax), SAVE :: a
  REAL(sp), SAVE :: epsold = -1.0_sp, xnew
  REAL(sp) :: eps1, errmax, fact, h, red, scale, wrkmin, xest
  REAL(sp), DIMENSION(SIZE(y)) :: yerr, ysav, yseq
  LOGICAL(lgt) :: reduct
  LOGICAL(lgt, SAVE :: first = .TRUE.

  ndum = assert_eq(SIZE(y), SIZE(dydx), SIZE(yscal), 'bsstep')
  IF (eps /= epsold) THEN              ! a new tolerance, reinitialize
     hnext = -1.0e29_sp                ! "impossible" values
     xnew=-1.0e29_sp
     eps1=safe1*eps
     a(:)=cumsum(nseq,1)
     WHERE (upper_triangle(kmaxx, kmaxx)) alf=eps1** &
          (outerdiff(a(2:), a(2:))/outerprod(arth( &
          3.0_sp, 2.0_sp, kmaxx), (a(2:)-a(1)+1.0_sp)))
     epsold=eps
     DO kopt=2,kmaxx-1                 ! determine optimal row number for
        IF (a(kopt+1) > a(kopt)*alf(kopt-1,kopt)) EXIT      ! convergence
     END DO
     kmax=kopt
  END IF
  h = htry
  ysav(:) = y(:)                       ! save the starting values
  IF (h /= hnext .OR. x /= xnew) THEN  ! a new stepsize or a new integration
     first = .TRUE.                    ! re-establish the order window
     kopt = kmax
  END IF
  reduct = .FALSE.
  main_loop: DO
     DO k = 1, kmax                    ! evaluate the sequence of modified
        xnew = x+h                     ! midpoint integrations
        IF (xnew == x) CALL nrerror('step size underflow in bsstep')
        CALL mmid(ysav, dydx, x, h, nseq(k), yseq, derivs)
        xest=(h/nseq(k))**2            !squared, since error series is even
        CALL pzextr(k, xest, yseq, y, yerr) ! perform extrapolation
        IF (k /= 1) THEN                  ! computer normalized error estimate
           errmax=MAXVAL(ABS(yerr(:)/yscal(:)))
           errmax=MAX(tiny, errmax)/eps   ! scale error relative to tolerance
           km=k-1
           err(km)=(errmax/safe1)**(1.0_sp/(2*km+1))
        END IF
        IF (k /= 1 .AND. (k >= kopt - 1 .OR. first)) THEN  ! in order window
           IF (errmax < 1.0) EXIT main_loop      ! converged
           IF(k == kmax .OR. k == kopt+1) THEN   ! check for possible step
              red = safe2/err(km)                ! size reduction
              EXIT
           ELSE IF (k == kopt) THEN
              IF (alf(kopt-1, kopt) < err(km)) THEN
                 red=1.0_sp/err(km)
                 EXIT
              END IF
           ELSE IF (kopt == kmax) THEN
              IF(alf(km, kmax-1) < err(km)) THEN
                 red = alf(km, kmax-1) * safe2/err(km)
                 EXIT
              END IF
           ELSE IF (alf(km, kopt) < err(km)) THEN
              red=alf(km,kopt=1) / err(km)
              EXIT
           END IF
        END IF
     END DO
     red=MAX(MIN(red, redmin), redmax)   ! reduce step size by at least
     h = h*red                           ! redmin and at most redmax
     reduct = .TRUE.
  END DO main_loop                       ! try again
  x = xnew                               ! successful step taken
  hdid=h
  first=.FALSE.
  kopt=1+iminloc(a(2:km+1)*MAX(err(1:km), scalmx))

! Compute optimal row for convergence and corresponding stepsize

  scale=MAX(err(kopt=1), scalmx)
  wrkmin=scale*a(kopt)
  hnext=h/scale
  IF(kopt >= k .AND. kopt /= kmax .AND. .NOT. reduct) THEN ! check for possible
     fact = MAX(scale/alf(kopt-1, kopt), scalmx)           ! order increase
     IF(a(kopt+1)*fact <=wrkmin) THEN                      ! but not if step
        hnext = h/fact                                     ! size just reduced
        kopt=kopt+1                                        
     END IF
  END IF

END SUBROUTINE bsstep