diff --git a/Solver/function.cpp b/Solver/function.cpp
index 1d288617f25f63a40059c8d14d6c9561948ff245..cd1e5465a42bc7f037616fb9f3e5aca4fb62c1e3 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -517,6 +517,33 @@ function *functionScaleNew(const function *f0, const double s) {
   return new functionScale (f0, s);
 }
 
+// functionMean
+
+class functionMean : public function {
+  fullMatrix<double> _f0;
+public:
+  functionMean(const function *f0) : function(f0->getNbCol()) {
+    setArgument (_f0, f0);
+  }
+  void call(dataCacheMap *m, fullMatrix<double> &val)
+  {
+    double mean;
+    for(int j = 0; j < val.size2(); j++) {
+      mean = 0;
+      for(int i = 0; i < val.size1(); i++)
+        mean += _f0(i, j);
+      mean /= (double) val.size1();
+      for(int i = 0; i < val.size1(); i++)
+        val(i, j) = mean;
+    }
+  }
+};
+
+function *functionMeanNew(const function *f0) {
+  return new functionMean (f0);
+}
+
+
 // functionCoordinates (get XYZ coordinates)
 
 class functionCoordinates : public function {
diff --git a/Solver/function.h b/Solver/function.h
index f9af049036971bfced069872d817f2863a996055..6a5ee284c1718bc39d774126c740030570882c0a 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -296,6 +296,7 @@ function *functionProdNew (const function *f0, const function *f1);
 function *functionScaleNew (const function *f0, const double s);
 function *functionExtractCompNew (const function *f0, const int iComp);
 function *functionCatCompNew(std::vector<const function *> fArray);
+function *functionMeanNew(const function *f0);
 
 function *getFunctionCoordinates();