anchor.diff:        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
  +        ! Increasing h to the minimum value,
           h = hmin
  -        ! and locally increasing the desired error (it will go back for the next step).
  +        ! and locally increasing the desired error, we prevent this from occurring.
  +        ! The error goes back to normal for the next integration step.
           der = der + .0001*der
         endif
         go to 1
  
positive.code:   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                  !
      !===========================================================!
      integer n, j
      double precision x,dy(n),yi(n),ys(n),yf(n),&
                       k1(n),k2(n),k3(n),k4(n),&
  
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),&
  