diff --git a/Solver/functionDerivator.cpp b/Solver/functionDerivator.cpp
index adbd0070bb7cbeddf48451fbd358e361f5cb2f31..4133c1bc712a668ceebb135e79a3b3f0442f4219 100644
--- a/Solver/functionDerivator.cpp
+++ b/Solver/functionDerivator.cpp
@@ -2,19 +2,20 @@
 #include "function.h"
 
 const fullMatrix<double> &functionDerivator::compute() {
-  _xRef = _x.get();
+  const int nbColX = _x.get().size2();
+  const double o_eps = 1. / _epsilon;
   _fRef = _f.get();
-  _dfdx.resize(_fRef.size1(),_fRef.size2()*_xRef.size2(), false);
-  for (int j=0;j<_xRef.size2();j++) {
+  _dfdx.resize(_fRef.size1(),_fRef.size2() * nbColX, false);
+  for (int j = 0; j < nbColX; ++j) {
     fullMatrix<double> &x =_x.set();
     for (int i=0;i<_fRef.size1();i++)
       x(i,j) += _epsilon;
     const fullMatrix<double> &f = _f.get();
     for (int k=0; k<_fRef.size2(); k++) 
       for (int i=0; i<_fRef.size1(); i++) 
-        _dfdx(i, k*_xRef.size2()+j) = (f(i,k)-_fRef(i,k))/_epsilon;
+        _dfdx(i, k * nbColX + j) = (f(i,k)-_fRef(i,k)) * o_eps;
     for (int i=0;i<_fRef.size1();i++)
-      x(i,j) = _xRef(i,j);
+      x(i,j) -= _epsilon;
   }
   _x.set();
   return _dfdx;