Updated NumericalDiff finite difference step so it is never smaller than epsilon

libeigen/eigen!2733

Closes #3105

unsupported/Eigen/src/NumericalDiff/NumericalDiff.h

diff --git a/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h b/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h
index a0b8c51..f07d3cc 100644
--- a/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h
+++ b/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h
@@ -87,10 +87,8 @@
 
     // Function Body
     for (int j = 0; j < n; ++j) {
-      h = eps * abs(x[j]);
-      if (h == 0.) {
-        h = eps;
-      }
+      const Scalar x_abs = abs(x[j]);
+      h = numext::maxi(x_abs, Scalar(1)) * eps;
       switch (mode) {
         case Forward:
           x[j] += h;

unsupported/test/NumericalDiff.cpp

diff --git a/unsupported/test/NumericalDiff.cpp b/unsupported/test/NumericalDiff.cpp
index b46dfb2..bfde6b6 100644
--- a/unsupported/test/NumericalDiff.cpp
+++ b/unsupported/test/NumericalDiff.cpp
@@ -97,7 +97,26 @@
   VERIFY_IS_APPROX(jac, actual_jac);
 }
 
+void test_small() {
+  VectorXd x(3);
+  MatrixXd jac(15, 3);
+  MatrixXd actual_jac(15, 3);
+  my_functor functor;
+
+  x << 0.082, 1.13, 1e-8;
+
+  // real one
+  functor.actual_df(x, actual_jac);
+
+  // using NumericalDiff
+  NumericalDiff<my_functor, Central> numDiff(functor);
+  numDiff.df(x, jac);
+
+  VERIFY_IS_APPROX(jac, actual_jac);
+}
+
 EIGEN_DECLARE_TEST(NumericalDiff) {
   CALL_SUBTEST(test_forward());
   CALL_SUBTEST(test_central());
+  CALL_SUBTEST(test_small());
 }