anchor.diff:    end function kdelta
   
     subroutine rk4gn(f,x,yi,yf,h,n)
  -    !===========================================================!
  -    ! Solve n first-order ODEs or n/2 second-order.             !
  -    ! Runge-Kutta 4 Gill's method                               !
  -    ! Daniel Celis Garza 24 Sept. 2015                          !
  -    !-----------------------------------------------------------!
  -    ! f(y,x,dx,n)   = ODEs to be solved                         !
  -    ! x             = independent variable @ step i             !
  -    ! yi()          = array of dependent variables @ step i     !
  -    ! ys()          = array of dependent variables @ stage s    !
  -    ! yf()          = array of dependent variables @ step i+1   !
  -    ! n             = ODE system dimension                      !
  -    !-----------------------------------------------------------!
  -    ! j             = counter variable                          !
  -    ! h             = step size                                 !
  -    ! ki()          = array of Runge-Kutta k's                  !
  -    !===========================================================!
  +    !================================================================!
  +    ! Solve n first-order ODEs or n/2 second-order.                  !
  +    ! Runge-Kutta 4 Gill's method                                    !
  +    ! Daniel Celis Garza 24 Sept. 2015                               !
  +    !----------------------------------------------------------------!
  +    ! f(x,y,dy,n) = ODEs to be solved                                !
  +    ! n           = ODE system dimension                             !
  +    ! h           = step size                                        !
  +    ! x           = independent variable @ step i                    !
  +    ! yi()        = array of dependent variables @ step i            !
  +    ! ys()        = array of dependent variables @ stage s of step i !
  +    ! yf()        = array of dependent variables @ step i+1          !
  +    ! er()        = array of integration errors                      !
  +    ! ki()        = array of Runge-Kutta k/h                         !
  +    ! j           = counter variable                                 !
  +    !================================================================!
       integer n, j
       double precision x,dy(n),yi(n),ys(n),yf(n),&
                        k1(n),k2(n),k3(n),k4(n),&
  
positive.code:       if (xtest == xi) then
          ! This means that the error is so strict, h is effectively zero,
          ! so the program will enter an infinite loop.
          ! This prevents it from doing so, by increasing h to the minimum value
          h = hmin
          ! and locally increasing the desired error (it will go back for the next step).
          der = der + .0001*der
        endif
        go to 1
  
negative.code:   end subroutine rk4gn
  
    subroutine rkckn(f,x,yi,yf,er,h,n)
      !================================================================!
      ! Solve n first-order ODEs or n/2 second-order.                  !
      ! Cash-Karp RK45                                                 !
      ! Daniel Celis Garza 24 Sept. 2015                               !
      !----------------------------------------------------------------!
      ! f(x,y,dy,n) = ODEs to be solved                                !
      ! n           = ODE system dimension                             !
      ! h           = step size                                        !
      ! x           = independent variable @ step i                    !
      ! yi()        = array of dependent variables @ step i            !
      ! ys()        = array of dependent variables @ stage s of step i !
      ! yf()        = array of dependent variables @ step i+1          !
      ! er()        = array of integration errors                      !
      ! ki()        = array of Runge-Kutta k/h                         !
      ! ci          = Butcher Table c-vector                           !
      ! aij         = Butcher Table A-matrix                           !
      ! bi          = Butcher Table b-vector                           !
      ! dbi         = b-b* vector difference for error calculation     !
      ! j           = counter variable                                 !
      !================================================================!
      integer n, j
      double precision x,dy(n),yi(n),ys(n),yf(n),&
                       k1(n),k2(n),k3(n),k4(n),k5(n),k6(n),er(n),&
  