diff --git a/Solver/function.cpp b/Solver/function.cpp
index c3dcb4dc4ef5753d925716c66bf52cdd787ee5e3..a2dda164a942db29ca7ed79443e3d0d92d4d2570 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -384,6 +384,40 @@ function *functionSumNew(const function *f0, const function *f1)
   return new functionSum (f0, f1);
 }
 
+// functionMinus
+
+class functionMinus : public function {
+ public:
+  fullMatrix<double> _f0, _f1;
+  void call(dataCacheMap *m, fullMatrix<double> &val) 
+  {
+    if (_f0.size2() != _f1.size2()) {
+      Msg::Error("trying to substract 2 functions of different sizes: %d - %d\n",
+                 _f0.size2(), _f1.size2());
+      throw;
+    }
+    for (int i = 0; i < val.size1(); i++)
+      for (int j = 0; j < val.size2(); j++)
+        val(i, j)= _f0(i, j) - _f1(i, j);
+  }
+  functionMinus(const function *f0, const function *f1) : function(f0->getNbCol()) 
+  {
+/*    if (f0->getNbCol() != f1->getNbCol()) {
+      Msg::Error("trying to substract 2 functions of different sizes: %d - %d\n",
+                 f0->getNbCol(), f1->getNbCol());
+      throw;
+    }*/
+    setArgument (_f0, f0);
+    setArgument (_f1, f1);
+  }
+};
+
+function *functionMinusNew(const function *f0, const function *f1) 
+{
+  return new functionMinus (f0, f1);
+}
+
+
 // functionLevelset
 
 class functionLevelset : public function {
diff --git a/Solver/function.h b/Solver/function.h
index 560af366d6fcc0348ca25a6e3138c2399a13dbfc..d8896a7c15bce85ca23831dc1f44bc70826da905 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -289,6 +289,7 @@ class functionC : public function {
 function *functionLevelsetNew (const function *f0, const double valMin, const double valPlus);
 function *functionLevelsetSmoothNew (const function *f0, const double valMin, const double valPlus, const double E);
 function *functionSumNew (const function *f0, const function *f1);
+function *functionMinusNew (const function *f0, const function *f1);
 function *functionProdNew (const function *f0, const function *f1);
 function *functionScaleNew (const function *f0, const double s);
 function *functionExtractCompNew (const function *f0, const int iComp);