diff --git a/Solver/TESTCASES/rect.geo b/Solver/TESTCASES/rect.geo
index 6f581ea1d64723210792af9058832d22aad94c69..f5124e89fc50c696cb1ea37b830004cb078d346f 100644
--- a/Solver/TESTCASES/rect.geo
+++ b/Solver/TESTCASES/rect.geo
@@ -21,5 +21,5 @@ Physical Line("Top") = {6};
 //Recombine Surface {9, 11};
 
 Field[1] = MathEval;
-Field[1].F = "0.01*1+0.01*1.7*y";
+Field[1].F = "0.01*1+0.01*100*(y*y)";
 Background Field = 1;
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index 2f63ec8fd30a3828ea68b9ba1dfd6bb40ff42a65..8abc615e793f4efeecaae1ba75d609c7c07d7c26 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -160,6 +160,13 @@ void dgDofContainer::scale(double f)
   _data->scale(f);
   _ghostData->scale(f); 
 }
+void dgDofContainer::scale(std::vector<dgGroupOfElements*>groupsVector, double f){
+  for(int i=0;i<groupsVector.size();i++){
+    dgGroupOfElements *g=groupsVector[i];
+    fullMatrix<double> &proxy=getGroupProxy(g);
+    proxy.scale(f);
+  }
+}
 
 void dgDofContainer::axpy(dgDofContainer &x, double a)
 {
diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h
index 6c3f15ac14f305b70a4ff8b3e8973921ce461a09..71395eda17fb4a2419b3652f6d30c0be1c7e8d81 100644
--- a/Solver/dgDofContainer.h
+++ b/Solver/dgDofContainer.h
@@ -26,6 +26,7 @@ private:
   int _mshStep;
 public:
   void scale(double f);
+  void scale(std::vector<dgGroupOfElements*>groups, double f);
   double norm();
   void axpy(dgDofContainer &x, double a=1.);
   void axpy(std::vector<dgGroupOfElements*>groups,dgDofContainer &x, double a=1.);
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 48b36559f1970c48d63340d4ebcbbb90fa32718d..837cd39486eb47031ac6d7ba8d27747520f70222 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -920,8 +920,9 @@ void dgGroupCollection::buildGroupsOfInterfaces() {
 
 
 // Split the groups of elements depending on their local time step
-void dgGroupCollection::splitGroupsForMultirate(double dtRef,dgConservationLaw *claw, dgDofContainer *solution){
+double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLaw *claw, dgDofContainer *solution){
   Msg::Info("Splitting Groups for multirate time stepping");
+  maxLevels--;
   int maxNumElems=getElementGroup(0)->getElement(0)->getGlobalNumber()+1;
   std::vector<int>oldGroupIds;
   oldGroupIds.resize(maxNumElems);
@@ -1040,19 +1041,19 @@ void dgGroupCollection::splitGroupsForMultirate(double dtRef,dgConservationLaw *
   dtMax=dtMax_max;
 #endif
   // dtMin is the time step for the most constrained element.
-  if(dtRef<=dtMin){
-    Msg::Info("No multirate, the reference time step is stable for all elements.");
-    return;
-  }
 
   Msg::Info("Time step for standard RK should be %e",dtMin);
   Msg::Info("Multirate base time step should be %e",dtMax);
 
-  dtRef=dtMax*0.8;
+  double dtRef=dtMax*0.8;
   // time steps are dtRef*2^(-dtExponent), with dtExponent ranging in [0:dtMaxExponent]
   int dtMaxExponent=0;
   while(dtRef/pow(2.0,(double)dtMaxExponent)>dtMin)
     dtMaxExponent++;
+  if(dtMaxExponent>maxLevels){
+    dtRef*=1.0/pow(2.0,(double)(dtMaxExponent-maxLevels));
+    dtMaxExponent=maxLevels;
+  }
   _dtMaxExponent=dtMaxExponent;
   std::vector<MElement *>currentNewGroup;
   std::vector< dgGroupOfElements* >newGroups;// indexed by newGroupId
@@ -1093,7 +1094,7 @@ void dgGroupCollection::splitGroupsForMultirate(double dtRef,dgConservationLaw *
         dgGroupOfElements *elGroup=getElementGroup(iGroup);
         for(int iElement=0;iElement<elGroup->getNbElements();iElement++){
           MElement *el=elGroup->getElement(iElement);
-          if(localDt[el->getNum()]>=currentDt && localDt[el->getNum()]<currentDt*2){
+          if(localDt[el->getNum()]>=currentDt && (localDt[el->getNum()]<currentDt*2 || currentExponent==0)){
             if(newGroupIds[el->getNum()]==-1){
               mapNewGroups[iGroup].push_back(el);
               newGroupIds[el->getNum()]=-2;
@@ -1194,6 +1195,7 @@ void dgGroupCollection::splitGroupsForMultirate(double dtRef,dgConservationLaw *
   Msg::Info("That makes a total of %d elements",count);
   _elementGroups.clear();
   _elementGroups=newGroups;
+  return dtRef;
 }
 
 
@@ -1243,7 +1245,7 @@ void dgGroupCollection::registerBindings(binding *b){
   cm->setDescription("Build the group of interfaces, i.e. boundary interfaces and inter-element interfaces");
   cm = cb->addMethod("splitGroupsForMultirate",&dgGroupCollection::splitGroupsForMultirate);
   cm->setDescription("Split the groups according to their own stable time step");
-  cm->setArgNames("refDt","claw","solution",NULL);
+  cm->setArgNames("maxLevels","claw","solution",NULL);
   cm = cb->addMethod("getNbElementGroups", &dgGroupCollection::getNbElementGroups);
   cm->setDescription("return the number of dgGroupOfElements");
   cm = cb->addMethod("getNbFaceGroups", &dgGroupCollection::getNbFaceGroups);
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index f108b15c06dbf89e40f098d002084e4ecff9844b..454403bbfca84816de6c25c0c0cd8d4982ad0c98 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -245,7 +245,7 @@ class dgGroupCollection {
   void buildGroupsOfElements (GModel *model,int dimension, int order);
   void buildGroupsOfInterfaces ();
 
-  void splitGroupsForMultirate(double refDt,dgConservationLaw *claw, dgDofContainer *solution);
+  double splitGroupsForMultirate(int maxLevels,dgConservationLaw *claw, dgDofContainer *solution);
 
   void find (MElement *elementToFind, int &iGroup, int &ithElementOfGroup);
 
diff --git a/Solver/dgRungeKutta.cpp b/Solver/dgRungeKutta.cpp
index a1b03cbdc3714b53888bd6adff1fcf517c0bfaea..683adc2695cbd41363180401e333df4b63d87aea 100644
--- a/Solver/dgRungeKutta.cpp
+++ b/Solver/dgRungeKutta.cpp
@@ -4,6 +4,7 @@
 #include "dgLimiter.h"
 #include "dgResidual.h"
 #include "dgGroupOfElements.h"
+#include <stdio.h>
 
 double dgRungeKutta::iterateEuler(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
   double A[] = {0};
@@ -423,7 +424,7 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
   _residualVolume=new dgResidualVolume(*_law);
   _residualInterface=new dgResidualInterface(*_law);
   _K=new dgDofContainer*[10];
-  for(int i=0;i>10;i++){
+  for(int i=0;i<10;i++){
     _K[i]=new dgDofContainer(gc,_law->getNbFields());
     _K[i]->setAll(0.0);
   }
@@ -476,6 +477,7 @@ dgRungeKuttaMultirate::~dgRungeKuttaMultirate(){
 }
 
 void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
+  double localDt=_dt/pow(2.0,(double)exponent);
   double _A[4][4]={
      {0, 0, 0, 0},
      {1.0/2.0, 0, 0 ,0},
@@ -494,8 +496,6 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
     {1.0/12.0,  1.0/6.0,   1.0/6.0,   1.0/12.0,  0,         1.0/6.0,   -1.0/6.0,  1.0/2.0,   0,         0},
     {1.0/12.0,  1.0/6.0,   1.0/6.0,   1.0/12.0,  0,         1.0/12.0,  1.0/6.0,   1.0/6.0,   1.0/12.0, 0}
   };
-  double bInner[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
-  double cInner[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
 
   // Big step RK43
   double _AOuter[10][10]={
@@ -510,29 +510,40 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
     {1.0/3.0 ,  0,         0,         0,         -1.0/3.0,  1,         0,         0,         0,         0},
     {1.0/3.0 ,  0,         0,         0,         -1.0/3.0,  1,         0,         0,         0,         0},
   };
-  double bOuter[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
-  double cOuter[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
   if(exponent>_maxExponent)
     return;
+//  Msg::Info("K%d, %d, %s",iK,exponent,isBuffer?"Buffer":"Bulk");
   // compute K[iK] for this exponent and buffer
-  _currentInput->scale(0.0);
   if(isBuffer){
-    std::vector<dgGroupOfElements *>gVi=_innerBufferGroupsOfElements[exponent].first;
-    std::vector<dgGroupOfElements *>gVo=_outerBufferGroupsOfElements[exponent].first;
+    std::vector<dgGroupOfElements *>&gVi=_innerBufferGroupsOfElements[exponent].first;
+    std::vector<dgGroupOfElements *>&gVo=_outerBufferGroupsOfElements[exponent].first;
+    _currentInput->scale(gVi,0.0);
+    _currentInput->scale(gVo,0.0);
     _currentInput->axpy(gVi,*_solution,1.0);
     _currentInput->axpy(gVo,*_solution,1.0);
     for(int i=0;i<iK;i++){
-      _currentInput->axpy(gVi,*_K[i],_AInner[iK][i]);
-      _currentInput->axpy(gVo,*_K[i],_AOuter[iK][i]);
+      if(_AInner[iK][i]!=0.0)
+        _currentInput->axpy(gVi,*_K[i],_AInner[iK][i]*localDt);
+      if(_AOuter[iK][i]!=0.0)
+        _currentInput->axpy(gVo,*_K[i],_AOuter[iK][i]*localDt);
     }
   }
   else{
-    std::vector<dgGroupOfElements *>gV=_bulkGroupsOfElements[exponent].first;
+    std::vector<dgGroupOfElements *> &gV=_bulkGroupsOfElements[exponent].first;
+    _currentInput->scale(gV,0.0);
     _currentInput->axpy(gV,*_solution,1.0);
     for(int i=0;i<iK;i++){
-      _currentInput->axpy(gV,*_K[i],_A[iK][i]);
+      if(_A[iK][i]!=0.0)
+        _currentInput->axpy(gV,*_K[i],_A[iK][i]*localDt);
     }
   }
+  computeResidual(exponent,isBuffer,_currentInput,_K[iK]);
+  /*
+  printf("%s%d\t",isBuffer?"B":"N",exponent);
+  for(int i=0;i<exponent;i++)
+    printf("    |");
+  printf(" K%d\n",iK);
+  */
   if(!isBuffer){
     switch(iK){
       case 0:
@@ -555,28 +566,55 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
   }
   else{
     if( (iK%5)<4 ){
-      computeK(iK%5,exponent+1,false);
+      computeK(iK%5,exponent,false);
     }
   }
   if( (iK==3 && !isBuffer) || (iK==9 && isBuffer) ){
     updateSolution(exponent,isBuffer);
   }
 }
-void dgRungeKuttaMultirate::updateSolution(int exponent,bool isBuffer){}
+void dgRungeKuttaMultirate::updateSolution(int exponent,bool isBuffer){
+  double localDt=_dt/pow(2.0,(double)exponent);
+  double b[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
+  double c[4]={0, 1.0/2.0, 1.0/2.0, 1};
+  double bInner[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
+  double cInner[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
+  double bOuter[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
+  double cOuter[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
+  if(isBuffer){
+    std::vector<dgGroupOfElements *>&gVi=_innerBufferGroupsOfElements[exponent].first;
+    std::vector<dgGroupOfElements *>&gVo=_outerBufferGroupsOfElements[exponent].first;
+    for(int i=0;i<10;i++){
+      if(bInner[i]!=0.0)
+        _solution->axpy(gVi,*_K[i],bInner[i]*localDt);
+      if(bOuter[i]!=0.0)
+        _solution->axpy(gVi,*_K[i],bOuter[i]*localDt);
+    }
+  }
+  else{
+    std::vector<dgGroupOfElements *>&gV=_bulkGroupsOfElements[exponent].first;
+    for(int i=0;i<4;i++){
+      if(b[i]!=0.0)
+        _solution->axpy(gV,*_K[i],b[i]*localDt);
+    }
+  }
+}
 void dgRungeKuttaMultirate::computeResidual(int exponent,bool isBuffer,dgDofContainer *solution, dgDofContainer *residual){
   if(isBuffer){
-    std::vector<dgGroupOfElements *>*vei=&_innerBufferGroupsOfElements[exponent].first;
-    std::vector<dgGroupOfElements *>*veo=&_outerBufferGroupsOfElements[exponent].first;
-    for(int i=0;i<vei->size();i++){
-      _residualVolume->computeAndMap1Group(*(*vei)[i], *solution, *residual);
+    std::vector<dgGroupOfElements *>&vei=_innerBufferGroupsOfElements[exponent].first;
+    std::vector<dgGroupOfElements *>&veo=_outerBufferGroupsOfElements[exponent].first;
+    residual->scale(vei,0.0);
+    residual->scale(veo,0.0);
+    for(int i=0;i<vei.size();i++){
+      _residualVolume->computeAndMap1Group(*vei[i], *solution, *residual);
     }
-    for(int i=0;i<veo->size();i++){
-      _residualVolume->computeAndMap1Group(*(*veo)[i], *solution, *residual);
+    for(int i=0;i<veo.size();i++){
+      _residualVolume->computeAndMap1Group(*veo[i], *solution, *residual);
     }
-    std::vector<dgGroupOfFaces *>* vfi=&_innerBufferGroupsOfElements[exponent].second;
-    std::vector<dgGroupOfFaces *>* vfo=&_outerBufferGroupsOfElements[exponent].second;
+    std::vector<dgGroupOfFaces *>& vfi=_innerBufferGroupsOfElements[exponent].second;
+    std::vector<dgGroupOfFaces *>& vfo=_outerBufferGroupsOfElements[exponent].second;
     std::set<dgGroupOfFaces *>gOFSet;
-    gOFSet.insert(vfo->begin(),vfo->end());
+    gOFSet.insert(vfo.begin(),vfo.end());
     for(std::set<dgGroupOfFaces *>::iterator it=gOFSet.begin();it!=gOFSet.end();it++){
       dgGroupOfFaces *faces=*it;
       const dgGroupOfElements *gL=&faces->getGroupLeft();
@@ -592,13 +630,14 @@ void dgRungeKuttaMultirate::computeResidual(int exponent,bool isBuffer,dgDofCont
     }
   }
   else{
-    std::vector<dgGroupOfElements *> *ve=&_bulkGroupsOfElements[exponent].first;
-    for(int i=0;i<ve->size();i++){
-      _residualVolume->computeAndMap1Group(*(*ve)[i], *solution, *residual);
+    std::vector<dgGroupOfElements *> &ve=_bulkGroupsOfElements[exponent].first;
+    residual->scale(ve,0.0);
+    for(int i=0;i<ve.size();i++){
+      _residualVolume->computeAndMap1Group(*ve[i], *solution, *residual);
     }
-    std::vector<dgGroupOfFaces *> *vf=&_bulkGroupsOfElements[exponent].second;
-    for(int i=0;i<vf->size();i++){
-      dgGroupOfFaces *faces=(*vf)[i];
+    std::vector<dgGroupOfFaces *> &vf=_bulkGroupsOfElements[exponent].second;
+    for(int i=0;i<vf.size();i++){
+      dgGroupOfFaces *faces=vf[i];
       const dgGroupOfElements *gL=&faces->getGroupLeft();
       const dgGroupOfElements *gR=&faces->getGroupRight();
       fullMatrix<double> solInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
@@ -615,6 +654,9 @@ void dgRungeKuttaMultirate::computeResidual(int exponent,bool isBuffer,dgDofCont
 
 double dgRungeKuttaMultirate::iterate(double dt, dgDofContainer *solution){
   _solution=solution;
+  _dt=dt;
+  _currentInput->scale(0.0);
+  _currentInput->axpy(*solution,1.0);
   computeK(0,0,false);
   computeK(1,0,false);
   computeK(2,0,false);
diff --git a/Solver/dgRungeKutta.h b/Solver/dgRungeKutta.h
index 4168dd6ad0fb61a6c0d4629969c25b5f4bfcc62e..fa17c931e4731353da8f829ebfc9a6ba54128ebb 100644
--- a/Solver/dgRungeKutta.h
+++ b/Solver/dgRungeKutta.h
@@ -34,6 +34,7 @@ class dgRungeKuttaMultirate: public dgRungeKutta{
   private:
   int _maxExponent;
   int _minExponent;
+  double _dt;
   dgConservationLaw *_law;
   dgDofContainer *_solution;
   dgDofContainer **_K;