anchor.diff:    integer :: m, n
     
     double precision, dimension(nspec), intent(in) :: Xin
  -  double precision, dimension(nspeC), intent(out) :: Xout
  +  double precision, dimension(nspec), intent(out) :: Xout
     double precision, intent(in) :: T, rho, tmax
     double precision, intent(out) :: enucdot
   
  
positive.code: 
       ! J_{m,n} = df_m/dX_n
       do n = 1, nspec
          X1(:) = Xin(:)
          X1(n) = X1(n) + eps
  
          X2(:) = Xin(:)
          X2(n) = Xin(n) - eps
  
          call rhs(rho, T, X1, dX1dt)
          call rhs(rho, T, X2, dX2dt)
  
negative.code:      dt = min(dt, Xin(n)/dXdt(n))
    enddo
  
    Xout(:) = xin(:)
    
    ! do the integration
    time = 0.0d0
    do while (time < tmax)
  