diff --git a/Solver/functionDerivator.cpp b/Solver/functionDerivator.cpp
index 79e99bccce080d8cb4855a01bc2f856ee904f92d..adbd0070bb7cbeddf48451fbd358e361f5cb2f31 100644
--- a/Solver/functionDerivator.cpp
+++ b/Solver/functionDerivator.cpp
@@ -1,19 +1,21 @@
 #include "functionDerivator.h"
 #include "function.h"
 
-void functionDerivator::compute() {
+const fullMatrix<double> &functionDerivator::compute() {
   _xRef = _x.get();
   _fRef = _f.get();
-  _dfdx = fullMatrix<double>(_fRef.size1(),_fRef.size2()*_xRef.size2());
+  _dfdx.resize(_fRef.size1(),_fRef.size2()*_xRef.size2(), false);
   for (int j=0;j<_xRef.size2();j++) {
-    _xDx = _xRef;
+    fullMatrix<double> &x =_x.set();
     for (int i=0;i<_fRef.size1();i++)
-      _xDx(i,j) += _epsilon;
-    _x.set()=_xDx;
+      x(i,j) += _epsilon;
     const fullMatrix<double> &f = _f.get();
-    for (int i=0; i<_fRef.size1(); i++) 
-      for (int k=0; k<_fRef.size2(); k++) 
+    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;
+    for (int i=0;i<_fRef.size1();i++)
+      x(i,j) = _xRef(i,j);
   }
-  _x.set()=_xRef;
+  _x.set();
+  return _dfdx;
 };
diff --git a/Solver/functionDerivator.h b/Solver/functionDerivator.h
index ad651e51604e168d849e071067f01cc9b7dee84c..1793a31930bea5936e2c00ead516200b3c32c737 100644
--- a/Solver/functionDerivator.h
+++ b/Solver/functionDerivator.h
@@ -6,7 +6,7 @@
 
 class functionDerivator {
   dataCacheDouble &_f,&_x;
-  fullMatrix<double> _fRef, _xRef,_xDx, _dfdx;
+  fullMatrix<double> _fRef, _xRef, _dfdx;
   double _epsilon;
   public:
   functionDerivator (dataCacheDouble &f, dataCacheDouble &x, double epsilon):
@@ -14,7 +14,7 @@ class functionDerivator {
     _x(x),
     _epsilon(epsilon)
     {}
-  void compute();
+  const fullMatrix<double> &compute();
   inline double get(int iQP, int iF, int iX) 
   {
     return _dfdx(iQP, iF*_xRef.size2()+iX);