diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 48a0b7a5dd469249a15e6213e317cfe1ccd23095..6a11a5bfadb414e48478a68506c0cddf328e5ac2 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -388,5 +388,10 @@ binding::binding()
 }
 
 binding *binding::_instance=NULL;
+binding::~binding() {
+  for (std::map<std::string,classBinding *>::iterator it =  classes.begin(); it != classes.end(); it++) {
+    delete it->second;
+  }
+}
 
 #endif
diff --git a/Common/LuaBindings.h b/Common/LuaBindings.h
index 8c95e3bde82f7db7bb58dbba2a6c6111faf851ed..f8a55bb9dfb570fec14e83e9b3f455b275650361 100644
--- a/Common/LuaBindings.h
+++ b/Common/LuaBindings.h
@@ -67,6 +67,7 @@ class binding {
   int readFile(const char *filename);
   void interactiveSession();
   binding();
+  ~binding();
   std::map<std::string,classBinding *> classes;
   template<class t>
   classBinding *addClass(std::string className);
@@ -1058,6 +1059,13 @@ class classBinding {
     return constructorLua;
   }
   inline const std::string getClassName()const {return _className;}
+  ~classBinding() {
+    if (_constructor)
+      delete _constructor;
+    for (std::map<std::string, luaMethodBinding *>::iterator it = methods.begin(); it!=methods.end(); it++) {
+      delete it->second;
+    }
+  }
 };
 
 template<typename t>
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 557c53c22d50cb06846bb7cfbb6f8ccdc4f3d9d6..2104c5bdde0c565e78cebeed15470c5c1223fc27 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -11,16 +11,19 @@
 #include "Bindings.h"
 
 function::~function() {
-  for (int i=0; i<arguments.size(); i++) {
-    delete arguments[i];
-  }
 }
+
 function::function(int nbCol, bool invalidatedOnElement):_nbCol(nbCol), _invalidatedOnElement(invalidatedOnElement){};
 
-functionReplace &function::addFunctionReplace() {
-  _functionReplaces.resize(_functionReplaces.size()+1);
-  return _functionReplaces.back();
+void function::addFunctionReplace(functionReplace &fr) {
+  _functionReplaces.push_back(&fr);
 }
+dataCacheDouble::~dataCacheDouble()
+{
+  for(int i = 0; i< functionReplaceCaches.size(); i++) {
+    delete functionReplaceCaches[i];
+  }
+};
 
 void dataCacheDouble::addMeAsDependencyOf (dataCacheDouble *newDep)
 {
@@ -59,57 +62,53 @@ dataCacheDouble::dataCacheDouble(dataCacheMap *m, function *f):
   }
   _dependencies.resize ( _function->arguments.size());
   for (unsigned int i=0;i<_function->arguments.size();i++) {
-    int iCache = _function->arguments[i]->iMap;
-    const function *f = _function->arguments[i]->f;
+    int iCache = _function->arguments[i].iMap;
+    const function *f = _function->arguments[i].f;
     _dependencies[i] = &m->getSecondaryCache(iCache)->get(f,this);
   }
   for (int i = 0; i < f->_functionReplaces.size(); i++) {
-    functionReplaceCaches.push_back (new functionReplaceCache(m, &f->_functionReplaces[i]));
+    functionReplaceCaches.push_back (new functionReplaceCache(m, f->_functionReplaces[i])); 
   }
 }
 
 void dataCacheDouble::resize() {
   _value = fullMatrix<double>(_nRowByPoint==0?1:_nRowByPoint*_cacheMap.getNbEvaluationPoints(),_value.size2());
 }
-
 void dataCacheDouble::_eval() {
-  for(unsigned int i=0;i<_substitutions.size(); i++){
+  /*for(unsigned int i=0;i<_substitutions.size(); i++){
     _substitutions[i].first->set() = (*_substitutions[i].second)();
-  }
+  }*/
   for(unsigned int i=0;i<_dependencies.size(); i++){
-    _function->arguments[i]->val.setAsProxy((*_dependencies[i])());
+    _function->arguments[i].val->setAsProxy((*_dependencies[i])());
   }
   for (int i = 0; i < _function->_functionReplaces.size(); i++) {
-    _function->_functionReplaces[i].currentCache = functionReplaceCaches[i];
+    _function->_functionReplaces[i]->currentCache = functionReplaceCaches[i];
     for (int j = 0; j < functionReplaceCaches[i]->toReplace.size() ; j++){
-      _function->_functionReplaces[i]._toReplace[j]->val.setAsProxy((*functionReplaceCaches[i]->toReplace[j])._value);
+      _function->_functionReplaces[i]->_toReplace[j].val->setAsProxy((*functionReplaceCaches[i]->toReplace[j])._value);
     }
   }
   _function->call(&_cacheMap, _value);
 }
 
 //dataCacheMap members
-dataCacheDouble &dataCacheMap::get(const function *f, dataCacheDouble *caller) 
-{
+dataCacheDouble &dataCacheMap::get(const function *f, dataCacheDouble *caller) {
   dataCacheDouble *&r = _cacheDoubleMap[f];
-  if (r==NULL && _parent) {
-    std::map<const function *, dataCacheDouble *>::iterator it = _parent->_cacheDoubleMap.find(f);
-    if (it != _parent->_cacheDoubleMap.end()) {
+  dataCacheMap *cParent = _parent;
+  while (cParent && r==NULL) {
+    std::map<const function *, dataCacheDouble *>::iterator it = cParent->_cacheDoubleMap.find(f);
+    if (it != cParent->_cacheDoubleMap.end()) {
       r = it->second;
       for (std::set<dataCacheDouble*>::iterator dep = r->_iDependOn.begin(); dep != r->_iDependOn.end(); dep++) {
         if (&(*dep)->_cacheMap == this) {
-          throw;
           r = NULL;
           break;
         }
       }
     }
+    cParent = cParent->_parent;
   }
-  if (r==NULL) {
-    if(_parent)
-      throw;
+  if (r==NULL)
     r = new dataCacheDouble(this, (function*)(f));
-  }
   if (caller)
     r->addMeAsDependencyOf(caller);
   return *r;
@@ -160,10 +159,38 @@ function *functionConstantNew(const std::vector<double> &v) {
   return new functionConstant(v);
 }
 
+
+
+class functionSum : public function {
+  public:
+  fullMatrix<double> _f0, _f1;
+  void call(dataCacheMap *m, fullMatrix<double> &val) {
+    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);
+      }
+  }
+  functionSum(const function *f0, const function *f1):function(f0->getNbCol()){
+    if (f0->getNbCol() != f1->getNbCol()) {
+      Msg::Error("trying to sum 2 functions of different sizes\n");
+      throw;
+    }
+    setArgument (_f0, f0);
+    setArgument (_f1, f1);
+  }
+};
+
+function *functionSumNew(const function *f0, const function *f1) {
+  return new functionSum (f0, f1);
+}
+
+
+
+
 // get XYZ coordinates
 class functionCoordinates : public function {
   static functionCoordinates *_instance;
-  const fullMatrix<double> &uvw;
+  fullMatrix<double> uvw;
   void call (dataCacheMap *m, fullMatrix<double> &xyz){
     for(int i = 0; i < uvw.size1(); i++){
       SPoint3 p;
@@ -173,9 +200,9 @@ class functionCoordinates : public function {
       xyz(i, 2) = p.z();
     }
   }
-  functionCoordinates():function(3),
-  uvw(addArgument(function::getParametricCoordinates()))
-    {
+  functionCoordinates():function(3)
+  {
+    setArgument(uvw, function::getParametricCoordinates());
   };// constructor is private only 1 instance can exists, call get to access the instance
  public:
   static function *get() {
@@ -249,7 +276,7 @@ function *function::getNormals() {
 }
 
 class functionStructuredGridFile : public function {
-  const fullMatrix<double> &coord;
+  fullMatrix<double> coord;
   public:
   int n[3];
   double d[3],o[3];
@@ -279,9 +306,8 @@ class functionStructuredGridFile : public function {
         +get(id[0]+1 ,id[1]+1 ,id[2]+1 )*(  xi[0])*(  xi[1])*(  xi[2]);
     }
   }
-  functionStructuredGridFile(const std::string filename, const function *coordFunction): function(1),
-     coord(addArgument(coordFunction))
-     {
+  functionStructuredGridFile(const std::string filename, const function *coordFunction): function(1) {
+    setArgument(coord, coordFunction);
     std::ifstream input(filename.c_str());
     if(!input)
       Msg::Error("cannot open file : %s",filename.c_str());
@@ -310,19 +336,21 @@ class functionStructuredGridFile : public function {
 class functionLua : public function {
   lua_State *_L;
   std::string _luaFunctionName;
+  std::vector<fullMatrix<double> > args;
  public:
   void call (dataCacheMap *m, fullMatrix<double> &res) {
     lua_getfield(_L, LUA_GLOBALSINDEX, _luaFunctionName.c_str());
     for (int i=0;i< arguments.size(); i++)
-      luaStack<const fullMatrix<double>*>::push(_L, &arguments[i]->val);
+      luaStack<const fullMatrix<double>*>::push(_L, &args[i]);
     luaStack<const fullMatrix<double>*>::push(_L, &res);
     lua_call(_L, arguments.size()+1, 0);
   }
   functionLua (int nbCol, std::string luaFunctionName, std::vector<const function*> dependencies, lua_State *L)
     : function(nbCol), _luaFunctionName(luaFunctionName), _L(L)
   {
-    for (std::vector<const function *>::iterator it = dependencies.begin(); it!= dependencies.end(); it++) {
-      addArgument(*it);
+    args.resize(dependencies.size());
+    for (int i = 0; i < dependencies.size(); i++) {
+      setArgument(args[i], dependencies[i]);
     }
   }
 };
@@ -343,50 +371,52 @@ void dataCacheMap::setNbEvaluationPoints(int nbEvaluationPoints) {
 
 //functionC
 class functionC : public function {
+  std::vector<fullMatrix<double> > args;
   void (*callback)(void);
   public:
   void call (dataCacheMap *m, fullMatrix<double> &val) {
-    switch (arguments.size()) {
+    switch (args.size()) {
       case 0 : 
         ((void (*)(fullMatrix<double> &))(callback))(val);
         break;
       case 1 : 
         ((void (*)(fullMatrix<double> &, const fullMatrix<double>&))
-         (callback)) (val, arguments[0]->val);
+         (callback)) (val, args[0]);
         break;
       case 2 : 
         ((void (*)(fullMatrix<double> &, const fullMatrix<double>&, const fullMatrix<double> &))
-         (callback)) (val, arguments[0]->val, arguments[1]->val);
+         (callback)) (val, args[0], args[1]);
         break;
       case 3 : 
         ((void (*)(fullMatrix<double> &, const fullMatrix<double>&, const fullMatrix<double>&, const fullMatrix<double>&))
-         (callback)) (val, arguments[0]->val, arguments[1]->val, arguments[2]->val);
+         (callback)) (val, args[0], args[1], args[2]);
         break;
       case 4 : 
         ((void (*)(fullMatrix<double> &, const fullMatrix<double>&, const fullMatrix<double>&, const fullMatrix<double>&,
                    const fullMatrix<double>&))
-         (callback)) (val, arguments[0]->val, arguments[1]->val, arguments[2]->val, arguments[3]->val);
+         (callback)) (val, args[0], args[1], args[2], args[3]);
         break;
       case 5 : 
         ((void (*)(fullMatrix<double> &, const fullMatrix<double>&, const fullMatrix<double>&, const fullMatrix<double>&,
                    const fullMatrix<double>&, const fullMatrix<double>&))
-         (callback)) (val, arguments[0]->val, arguments[1]->val, arguments[2]->val, arguments[3]->val, arguments[4]->val);
+         (callback)) (val, args[0], args[1], args[2], args[3], args[4]);
         break;
       case 6 : 
         ((void (*)(fullMatrix<double> &, const fullMatrix<double>&, const fullMatrix<double>&, const fullMatrix<double>&,
                    const fullMatrix<double>&, const fullMatrix<double>&, const fullMatrix<double>&))
-         (callback)) (val, arguments[0]->val, arguments[1]->val, arguments[2]->val, arguments[3]->val, arguments[4]->val, arguments[5]->val);
+         (callback)) (val, args[0], args[1], args[2], args[3], args[4], args[5]);
         break;
       default :
-        Msg::Error("C callback not implemented for %i argurments", arguments.size());
+        Msg::Error("C callback not implemented for %i argurments", args.size());
     }
   }
   functionC (std::string file, std::string symbol, int nbCol, std::vector<const function *> dependencies):
     function(nbCol)
   {
 #if defined(HAVE_DLOPEN)
-    for (std::vector<const function *>::iterator it = dependencies.begin(); it!= dependencies.end(); it++) {
-      addArgument(*it);
+    args.resize(dependencies.size());
+    for(int i=0; i < dependencies.size(); i++) {
+      setArgument(args[i], dependencies[i]);
     }
     void *dlHandler;
     dlHandler = dlopen(file.c_str(),RTLD_NOW);
@@ -461,35 +491,36 @@ void function::registerBindings(binding *b){
 #endif
 }
 
-fullMatrix<double> &functionReplace::replace(const function *f, int iMap) {
-  function::argument *arg = new function::argument(iMap, f);
-  _toReplace.push_back(arg);
-  return arg->val;
+void functionReplace::replace(fullMatrix<double> &v, const function *f, int iMap) {
+  _toReplace.push_back(function::argument(v, iMap, f));
 }
 
-const fullMatrix<double> &functionReplace::get(const function *f, int iMap) {
-  function::argument *arg = new function::argument(iMap, f);
-  _toCompute.push_back(arg);
-  return arg->val;
+void functionReplace::get(fullMatrix<double> &v, const function *f, int iMap) {
+  _toCompute.push_back(function::argument(v, iMap, f));
 }
 
 void functionReplace::compute(){
-  for (int i = 0; i < _toReplace.size(); i++) {
-  //printf("a %p\n",&currentCache->toReplace[i]->_value);
-    currentCache->toReplace[i]->set();//.setAsProxy(_toReplace[i]->val);
-  //printf("b\n");
-  }
+  for (int i = 0; i < _toReplace.size(); i++)
+    currentCache->toReplace[i]->set();
   for (int i = 0; i < _toCompute.size(); i++) {
-    _toCompute[i]->val.setAsProxy((*currentCache->toCompute[i])());
+    _toCompute[i].val->setAsProxy((*currentCache->toCompute[i])());
   }
 };
 
+
 functionReplaceCache::functionReplaceCache(dataCacheMap *m, functionReplace *rep) {
   map = m->newChild();
+  for (int i = 0; i < m->_secondaryCaches.size(); i ++) {
+    map->addSecondaryCache (m->getSecondaryCache(i+1)->newChild());
+  }
   for (int i = 0; i < rep->_toReplace.size(); i++) {
-    toReplace.push_back (&map->getSecondaryCache(rep->_toReplace[i]->iMap)->substitute(rep->_toReplace[i]->f));
+    toReplace.push_back (&map->getSecondaryCache(rep->_toReplace[i].iMap)->substitute(rep->_toReplace[i].f));
   }
   for (int i = 0; i < rep->_toCompute.size(); i++) {
-    toCompute.push_back (&map->getSecondaryCache(rep->_toCompute[i]->iMap)->get(rep->_toCompute[i]->f));
+    dataCacheMap *m2 = map->getSecondaryCache(rep->_toCompute[i].iMap);
+    toCompute.push_back (&m2->get(rep->_toCompute[i].f));
   }
 }
+functionReplaceCache::~functionReplaceCache() {
+  delete map;
+}
diff --git a/Solver/function.h b/Solver/function.h
index dd36dec5ce3916902325205b5bfc28ef6e79ef41..a7242f4fa2173794c8fe76a250e21ca94df060e3 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -26,13 +26,14 @@ class function {
     public:
     int iMap;
     const function *f;
-    fullMatrix<double> val;
-    argument(int iMap_, const function *f_){
+    fullMatrix<double> *val;
+    argument(fullMatrix<double> &v, int iMap_, const function *f_){
+      val = &v;
       iMap = iMap_;
       f = f_;
     }
   };
-  std::vector<functionReplace> _functionReplaces;
+  std::vector<functionReplace*> _functionReplaces;
   class substitutedFunction {
     public:
     int iMap;
@@ -43,12 +44,11 @@ class function {
   std::vector<int> _childrenCache;
   std::vector<substitutedFunction> _substitutedFunctions;
   virtual void call (dataCacheMap *m, fullMatrix<double> &res)=0;
-  std::vector<argument*> arguments;
-  const fullMatrix<double> &addArgument(const function *f, int iMap = 0) {
+  std::vector<argument> arguments;
+  const void setArgument(fullMatrix<double> &v, const function *f, int iMap = 0) {
     if(f==NULL)
       throw;
-    arguments.push_back(new argument(iMap, f));
-    return arguments.back()->val;
+    arguments.push_back(argument(v, iMap, f));
   }
   void addChildDataCacheMap(int parent) {
     _childrenCache.push_back(parent);
@@ -65,7 +65,7 @@ class function {
   function(int nbCol, bool invalidatedOnElement = true);
   inline int getNbCol()const {return _nbCol;}
   inline bool isInvalitedOnElement() { return _invalidatedOnElement;}
-  functionReplace &addFunctionReplace();
+  void addFunctionReplace(functionReplace &fr);
   
   static function *getSolution();
   static function *getSolutionGradient();
@@ -108,7 +108,6 @@ public :
 
   int _nRowByPoint;
   function *_function;
- protected:
   dataCacheMap &_cacheMap;
   // do the actual computation and put the result into _value
   // still virtual because it is overrided by conservation law terms, as soon as conservation law terms will be regular functions, we will remove this
@@ -149,7 +148,7 @@ public :
   void resize();
   dataCacheDouble(dataCacheMap *,function *f);
   dataCacheDouble(dataCacheMap &m, int nRowByPoint, int nbCol);
-  virtual ~dataCacheDouble(){};
+  virtual ~dataCacheDouble();
 };
 
 
@@ -188,9 +187,6 @@ class dataCacheMap {
   }
   void addSecondaryCache(dataCacheMap *s) {
     _secondaryCaches.push_back(s);
-    if(_secondaryCaches.size()>1){
-      printf("!!!!!!!!!!!!!!!!!!!!!\n");
-    }
   }
   dataCacheDouble &get(const function *f, dataCacheDouble *caller=0);
   dataCacheDouble &substitute(const function *f);
@@ -226,11 +222,11 @@ class functionReplace {
   friend class dataCacheDouble;
   protected:
   functionReplaceCache *currentCache;
-  std::vector <function::argument*> _toReplace;
-  std::vector <function::argument*> _toCompute;
+  std::vector <function::argument> _toReplace;
+  std::vector <function::argument> _toCompute;
   public :
-  const fullMatrix<double> &get(const function *, int iMap = 0);
-  fullMatrix<double> &replace(const function *, int iMap = 0);
+  void get(fullMatrix<double> &v, const function *, int iMap = 0);
+  void replace(fullMatrix<double> &v, const function *, int iMap = 0);
   void compute ();
 };
 
@@ -240,11 +236,13 @@ class functionReplaceCache {
   std::vector <dataCacheDouble*> toReplace;
   std::vector <dataCacheDouble*> toCompute;
   functionReplaceCache(dataCacheMap *m, functionReplace *rep);
+  ~functionReplaceCache();
 };
 
 
 function *functionConstantNew(const std::vector<double>&);
 function *functionConstantNew(double);
+function *functionSumNew (const function *f0, const function *f1);
 
 class functionSolution : public function {
   static functionSolution *_instance;