diff --git a/Apps/Common/SIMCoupledSI.h b/Apps/Common/SIMCoupledSI.h
index 36d65b11..fcb24f27 100644
--- a/Apps/Common/SIMCoupledSI.h
+++ b/Apps/Common/SIMCoupledSI.h
@@ -14,12 +14,14 @@
 #ifndef SIM_COUPLED_SI_H_
 #define SIM_COUPLED_SI_H_
 
-#include "IFEM.h"
-#include "Function.h"
-#include "MatVec.h"
 #include "SIMCoupled.h"
 #include "SIMenums.h"
+#include "MatVec.h"
 #include "TimeStep.h"
+#include "Functions.h"
+#include "Utilities.h"
+#include "IFEM.h"
+#include <tinyxml.h>
 
 
 /*!
@@ -33,20 +35,17 @@ public:
   //! \brief The constructor forwards to the parent class constructor.
   SIMCoupledSI(T1& s1, T2& s2) : SIMCoupled<T1,T2>(s1,s2)
   {
-    maxIter = -1;
-    maxIter0 = 50;
-    aitken = false;
+    maxIter = 50;
+    maxItFunc = nullptr;
+    noStg = aitken = false;
     omega = omega0 = 0.0;
   }
 
-  //! \brief Empty destructor.
-  virtual ~SIMCoupledSI() {}
+  //! \brief The destructor deletes the max iteration function.
+  virtual ~SIMCoupledSI() { delete maxItFunc; }
 
   //! \brief Enable/disable the staggering iteration cycles.
-  virtual void enableStaggering(bool enable = true)
-  {
-    maxIter = enable ? maxIter0 : 0;
-  }
+  virtual void enableStaggering(bool enable = true) { noStg = !enable; }
 
   //! \brief Returns residual to use for aitken acceleration.
   virtual const Vector& getAitkenResidual() const
@@ -62,7 +61,7 @@ public:
     return empty;
   }
 
-  //! \brief Set the relaxed solution.
+  //! \brief Sets the relaxed solution.
   virtual void setRelaxedSolution(const Vector&) {}
 
   //! \brief Computes the solution for the current time step.
@@ -72,11 +71,9 @@ public:
       this->S1.getProcessAdm().cout <<"\n  step="<< tp.step
                                     <<"  time="<< tp.time.t << std::endl;
 
-    if (maxIter == -1)
-      maxIter = maxIter0;
-
+    int maxit = this->getMaxit(tp.step);
     SIM::ConvStatus conv = SIM::OK;
-    for (tp.iter = 0; tp.iter <= this->getMaxit(tp.step) && conv != SIM::CONVERGED; tp.iter++)
+    for (tp.iter = 0; tp.iter <= maxit && conv != SIM::CONVERGED; tp.iter++)
     {
       SIM::ConvStatus status1 = SIM::OK, status2 = SIM::OK;
       if (firstS1 && (status1 = this->S1.solveIteration(tp)) <= SIM::DIVERGED)
@@ -91,24 +88,22 @@ public:
       if ((conv = this->checkConvergence(tp,status1,status2)) <= SIM::DIVERGED)
         return false;
 
-      // Calculate aitken acceleration factor
-      if (aitken && omega0 != 0.0) {
+      if (omega0 != 0.0) {
         if (tp.iter > 0) {
+          // Calculation of Aitken acceleration factor
+          if (aitken) {
             Vector r1 = this->getAitkenResidual();
             r1 -= prevRes;
             omega *= -prevRes.dot(r1) / r1.dot(r1);
             if (fabs(omega) < 1e-6) {
-              std::cerr << "\n** relaxation weight too small, resetting to default";
+              std::cerr <<"\n  ** Relaxation weight "<< omega <<" too small,"
+                        <<" resetting to default "<< omega0 << std::endl;
               omega = omega0;
             }
-        }
-        prevRes = this->getAitkenResidual();
-      }
+          }
 
-      // Perform relaxation
-      if (omega0 != 0.0) {
-        if (tp.iter > 0) {
-          IFEM::cout << ", omega=" << omega;
+          // Perform relaxation
+          IFEM::cout <<", omega="<< omega;
           prevSol *= 1.0 - omega;
           prevSol.add(this->getRelaxationVector(), omega);
           this->setRelaxedSolution(prevSol);
@@ -116,6 +111,8 @@ public:
           omega = omega0;
           prevSol = this->getRelaxationVector();
         }
+        if (aitken)
+          prevRes = this->getAitkenResidual();
       }
       IFEM::cout << std::endl;
     }
@@ -145,24 +142,53 @@ public:
   }
 
 protected:
-  //! \brief Returns the maximum number of iterations.
+  //! \brief Parses sub-iteration setup from an XML tag.
+  void parseIterations(const TiXmlElement* elem)
+  {
+    IFEM::cout <<"\tUsing sub-iterations\n";
+
+    std::string func;
+    if (utl::getAttribute(elem,"maxFunc",func) ||
+        utl::getAttribute(elem,"max",func))
+      if (func.find_first_of('t') != std::string::npos)
+      {
+        IFEM::cout <<"\t\tmax = ";
+        maxItFunc = utl::parseTimeFunc(func.c_str(),"expression");
+      }
+
+    if (!maxItFunc && utl::getAttribute(elem,"max",maxIter))
+      IFEM::cout <<"\t\tmax = "<< maxIter << std::endl;
+
+    if ((utl::getAttribute(elem,"relax",omega) ||
+         utl::getAttribute(elem,"omega",omega)) && omega != 0.0) {
+      IFEM::cout <<"\t\trelaxation = "<< omega;
+      if (utl::getAttribute(elem,"aitken",aitken) && aitken)
+        IFEM::cout <<" (aitken)";
+      IFEM::cout << std::endl;
+      omega0 = omega;
+    }
+  }
+
+  //! \brief Returns the maximum number of sub-iteration cycles.
   int getMaxit(int iStep = 0) const
   {
-    if (maxSubItFunc) {
-      Vec4 X;
-      X.t = iStep;
-      return static_cast<int>((*maxSubItFunc)(X));
-    }
+    if (noStg)
+      return 0; // staggering is disabled
+    else if (maxItFunc)
+      return static_cast<int>((*maxItFunc)(iStep));
 
     return maxIter;
   }
 
-  int maxIter0; //!< Initial maximum number of iterations
+private:
   int maxIter; //!< Maximum number of iterations
-  std::unique_ptr<RealFunc> maxSubItFunc; //!< Maximum number of sub-iterations as a function
-  double omega; //!< Relaxation parameter
+  ScalarFunc* maxItFunc; //!< Maximum number of iterations as a function
+  bool noStg;  //!< If \e true, sub-iterations is disabled
+
   double omega0; //!< Initial relaxation parameter
-  bool aitken; //!< True to enable aitken-acceleration
+  double omega;  //!< Current relaxation parameter
+  bool   aitken; //!< If \e true, Aitken-acceleration is enabled
+
   Vector prevSol; //!< Previous solution for relaxed field
   Vector prevRes; //!< Previous residual used for aitken-acceleration factor
 };