anchor.diff:       dt = min(dt, Xin(n)/dXdt(n))
     enddo
   
  +  Xout(:) = xin(:)
  +  
     ! do the integration
     time = 0.0d0
     do while (time < tmax)
  
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:   
    double precision :: dXdt(nspec)
    double precision :: X1(nspec), X2(nspec), dX1dt(nspec), dX2dt(nspec)
    double precision :: A(nspec, nspec), J(nspec, nspec), b(nspec)
    double precision :: dfdX
    integer :: info
    
    double precision, parameter :: eps = 1.d-8
  
    integer :: ipiv(nspec)
  
    ! get an estimate of the timestep by looking at the RHS
    ! dt = min{X/(dX/dt)}
  