anchor.diff:    real(kind=dp_t) :: time, dt, I
   
     integer :: m, n
  -    
  +
     real(kind=dp_t) :: dXdt(nspec)
     real(kind=dp_t) :: X1(nspec), X2(nspec), dX1dt(nspec), dX2dt(nspec)
     real(kind=dp_t) :: X_n(nspec), X_np1(nspec), dX(nspec), dfdX(nspec)
     real(kind=dp_t) :: A(nspec, nspec), J(nspec, nspec), b(nspec)
     integer :: info
  -  
  +
     real(kind=dp_t) :: eps = 1.d-8
   
     real(kind=dp_t) :: rpar(n_rpar_comps)
  
positive.code: 
          ! initial guess
          X_np1(:) = X_n(:)
       
          do iter = 1, max_iter
  
             ! construct a numerical Jacobian via simple differencing
             ! Lapack matrix ordering is described here:
             ! http://www.netlib.org/lapack/lug/node116.html#subsecarrayargs
             ! in general, A(i,j) holds matrix element a_{i,j}
          
             ! J_{m,n} = df_m/dX_n
             do n = 1, nspec
                X1(:) = X_np1(:)
                X1(n) = X1(n) + eps
                
                X2(:) = X_np1(:)
                X2(n) = X2(n) - eps
                
                call rhs(nspec, ZERO, X1, dX1dt, rpar, ipar)
                call rhs(nspec, ZERO, X2, dX2dt, rpar, ipar)
  
  
negative.code:                  else
                      I = 0.0d0
                   endif
  
                   A(m,n) = I - dt*J(m,n)
                enddo
             enddo
  
             ! construct the RHS
             call rhs(nspec, ZERO, X_np1, dXdt, rpar, ipar)
             b(:) = X_n(:) - X_np1(:) + dt*dXdt(:)
  
             ! solve the linear system -- we'll just use Lapack.  Note that the
             ! input b is overwritten on output
             call dgesv(nspec, 1, A, nspec, ipiv, b, nspec, info)
  
             dX(:) = b(:)
  
             X_np1(:) = X_np1(:) + dX(:)
  
             ! compare to see if we converged
  