anchor.diff:    end subroutine rk4gn
   
     subroutine rkckn(f,x,yi,yf,er,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.                  !
  +    ! 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),&
  
positive.code:   end subroutine rkckn
  
    subroutine rkcka(f,xi,xf,yi,yf,der,h,hmin,n)
      external f
      integer n, j
      double precision xi,xf,xtest,yi(n),yf(n),dy(n),yscal(n),der,er(n),h,ht,hn,hmin,maxer,miner,&
                       tiny,s,shrink,grow,ercor
      parameter (tiny=1e-30,s=0.9,shrink=-0.25,grow=-0.2,ercor=1.89e-4)
      call f(xi, yi, dy, n)
  
negative.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.
          ! Increasing h to the minimum value,
          h = hmin
          ! 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
  