diff --git a/unsupported/Eigen/src/NonLinear/MathFunctions.h b/unsupported/Eigen/src/NonLinear/MathFunctions.h
index f4688cc7d..4fedfdc9f 100644
--- a/unsupported/Eigen/src/NonLinear/MathFunctions.h
+++ b/unsupported/Eigen/src/NonLinear/MathFunctions.h
@@ -27,4 +27,22 @@
 
 #include <cminpack.h>
 
+template<typename Functor, typename VectorType>
+// TODO : fixe Scalar here
+int ei_hybrd1(
+        VectorType &x,
+        VectorType &fvec,
+//        ei_traits<VectorType>::Scalar tol
+        double  tol
+//        = ei::sqrt(machine_epsilon<VectorType::Scalar>())
+        )
+{
+    typedef typename VectorType::Scalar Scalar;
+    int lwa = (x.size()*(3*x.size()+13))/2;
+    VectorType wa(lwa);
+    fvec.resize(x.size());
+    return hybrd1(Functor::f, 0, x.size(), x.data(), fvec.data(), tol, wa.data(), lwa);
+}
+
+
 #endif // EIGEN_NONLINEAR_MATHFUNCTIONS_H
diff --git a/unsupported/test/NonLinear.cpp b/unsupported/test/NonLinear.cpp
index 3f8562841..105199fba 100644
--- a/unsupported/test/NonLinear.cpp
+++ b/unsupported/test/NonLinear.cpp
@@ -483,32 +483,36 @@ int fcn_hybrd1(void * /*p*/, int n, const double *x, double *fvec, int /*iflag*/
   return 0;
 }
 
+
+struct myfunctor {
+    static int f(void *p, int n, const double *x, double *fvec, int iflag )
+    { return fcn_hybrd1(p,n,x,fvec,iflag) ; }
+};
+
 void testHybrd1()
 {
-  int j, n, info, lwa;
-  double tol, fnorm;
-  double x[9], fvec[9], wa[180];
+  int j, n=9, info;
+  double fnorm;
+  Eigen::VectorXd x(9), fvec(9);
 
-  n = 9;
 
 /*      the following starting values provide a rough solution. */
 
-  for (j=1; j<=9; j++)
+  for (j=1; j<=n; j++)
     {
       x[j-1] = -1.;
     }
 
-  lwa = 180;
-
 /*      set tol to the square root of the machine precision. */
 /*      unless high solutions are required, */
 /*      this is the recommended setting. */
 
-  tol = sqrt(dpmpar(1));
-  info = hybrd1(fcn_hybrd1, 0, n, x, fvec, tol, wa, lwa);
-  fnorm = enorm(n, fvec);
+  info = ei_hybrd1<myfunctor,VectorXd>(x, fvec, sqrt(dpmpar(1)));
+
+  fnorm = enorm(fvec.size(), fvec.data());
+
+  VERIFY_IS_APPROX(fvec.norm(), 1.192636e-08);
 
-  VERIFY_IS_APPROX(fnorm, 1.192636e-08);
   VERIFY(info==1);
 
   double x_ref[] = {