diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/epsilon.h b/OMCompiler/SimulationRuntime/c/simulation/solver/epsilon.h
index 6cfea4ccea2..0470b959e9a 100644
--- a/OMCompiler/SimulationRuntime/c/simulation/solver/epsilon.h
+++ b/OMCompiler/SimulationRuntime/c/simulation/solver/epsilon.h
@@ -52,6 +52,7 @@ static const double DASSL_STEP_EPS = 1e-13;
  * defines the minimal step size
  */
 static const double MINIMAL_STEP_SIZE = 1e-12;
+static const double GB_MINIMAL_STEP_SIZE = 1e-20;
 
 /*
  * used in model_help.c for function setZCtol
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_conf.c b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_conf.c
index 84045a76652..d10db0a32c5 100644
--- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_conf.c
+++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_conf.c
@@ -58,6 +58,7 @@ const char *GB_INTERPOL_METHOD_NAME[GB_INTERPOL_MAX] = {
   /* GB_INTERPOL_UNKNOWN */           "unknown",
   /* GB_INTERPOL_LIN */               "linear",
   /* GB_INTERPOL_HERMITE */           "hermite",
+  /* GB_INTERPOL_HERMITE_a */         "hermite_a",
   /* GB_INTERPOL_HERMITE_b */         "hermite_b",
   /* GB_INTERPOL_HERMITE_ERRCTRL */   "hermite_errctrl",
   /* GB_DENSE_OUTPUT */               "dense_output",
@@ -68,7 +69,8 @@ const char *GB_INTERPOL_METHOD_DESC[GB_INTERPOL_MAX] = {
   /* GB_INTERPOL_UNKNOWN */         "unknown",
   /* GB_INTERPOL_LIN */             "Linear interpolation (1st order)",
   /* GB_INTERPOL_HERMITE */         "Hermite interpolation (3rd order)",
-  /* GB_INTERPOL_HERMITE_b */       "Hermite interpolation (only for left hand side)",
+  /* GB_INTERPOL_HERMITE_a */       "Hermite interpolation (only for left hand side)",
+  /* GB_INTERPOL_HERMITE_b */       "Hermite interpolation (only for right hand side)",
   /* GB_INTERPOL_HERMITE_ERRCTRL */ "Hermite interpolation with error control",
   /* GB_DENSE_OUTPUT */             "use dense output formula for interpolation",
   /* GB_DENSE_OUTPUT_ERRCTRL */     "use dense output fomular with error control"
@@ -77,9 +79,9 @@ const char *GB_INTERPOL_METHOD_DESC[GB_INTERPOL_MAX] = {
 /**
  * @brief Get Runge-Kutta method from simulation flag FLAG_SR or FLAG_MR.
  *
- * Defaults to method RK_LOBA_IIIB_4 for single-rate.
+ * Defaults to method RK_ESDIRK4 for single-rate.
  *
- * Defaults to method RK_SDIRK2 for multi-rate method, if single-rate method is implicit.
+ * Defaults to method RK_ESDIRK4 for multi-rate method, if single-rate method is implicit.
  * Otherwise us same method as single-rate method.
  *
  * Returns GB_UNKNOWN if flag is not known.
@@ -131,15 +133,15 @@ enum GB_METHOD getGB_method(enum _FLAG flag) {
     case RK_LOBA_IIIC_4:
       // Default value for inner integration method
       // if the outer integration method is full implicit
-      return RK_ESDIRK3;
+      return RK_ESDIRK4;
     default:
       return singleRateMethod;
     }
   }
 
   // Default value for single-rate method
-  infoStreamPrint(LOG_SOLVER, 0, "Chosen gbode method: esdirk3 [default]");
-  return RK_ESDIRK3;
+  infoStreamPrint(LOG_SOLVER, 0, "Chosen gbode method: esdirk4 [default]");
+  return RK_ESDIRK4;
 }
 
 /**
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_conf.h b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_conf.h
index 362503891d6..61ae661d4ac 100644
--- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_conf.h
+++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_conf.h
@@ -56,14 +56,14 @@ extern const char *GB_CTRL_METHOD_NAME[GB_CTRL_MAX];
 extern const char *GB_CTRL_METHOD_DESC[GB_CTRL_MAX];
 
 enum GB_INTERPOL_METHOD {
-  GB_INTERPOL_UNKNOWN = 0,          /* Unknown interpolation method */
-  GB_INTERPOL_LIN = 1,              /* Linear interpolation */
-  GB_INTERPOL_HERMITE = 2,          /* Hermite interpolation */
-  GB_INTERPOL_HERMITE_b = 3,        /* Hermite interpolation (only for left hand side)*/
-  GB_INTERPOL_HERMITE_ERRCTRL = 4,  /* Hermite interpolation with error control */
-  GB_DENSE_OUTPUT = 5,              /* Dense output, if available else hermite */
-  GB_DENSE_OUTPUT_ERRCTRL = 6,       /* Dense output, if available else hermite with error control */
-
+  GB_INTERPOL_UNKNOWN = 0,      /* Unknown interpolation method */
+  GB_INTERPOL_LIN,              /* Linear interpolation */
+  GB_INTERPOL_HERMITE,          /* Hermite interpolation */
+  GB_INTERPOL_HERMITE_a,        /* Hermite interpolation (only for left hand side)*/
+  GB_INTERPOL_HERMITE_b,        /* Hermite interpolation (only for right hand side)*/
+  GB_INTERPOL_HERMITE_ERRCTRL,  /* Hermite interpolation with error control */
+  GB_DENSE_OUTPUT,              /* Dense output, if available else hermite */
+  GB_DENSE_OUTPUT_ERRCTRL,      /* Dense output, if available else hermite with error control */
 
   GB_INTERPOL_MAX
 };
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_ctrl.c b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_ctrl.c
index 5d16ae0a856..abe629e361c 100644
--- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_ctrl.c
+++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_ctrl.c
@@ -227,7 +227,7 @@ void getInitStepSize(DATA* data, threadData_t* threadData, DATA_GBODE* gbData)
     h1 = fmax(1e-6, h0*1e-3);
   }
 
-  gbData->stepSize = 0.5*fmin(100*h0,h1)*50;
+  gbData->stepSize = 0.5*fmin(100*h0,h1);
   gbData->lastStepSize = 0.0;
 
   sData->timeValue = gbData->time;
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.c b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.c
index ca0d63f9705..2fb0c822738 100644
--- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.c
+++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_main.c
@@ -96,7 +96,7 @@ void gbode_fODE(DATA *data, threadData_t *threadData, unsigned int* counter)
  * @param solverInfo    Information about main solver.
  * @return int          Return 0 on success, -1 on failure.
  */
-int gbodef_allocateData(DATA *data, threadData_t *threadData, DATA_GBODE *gbData)
+int gbodef_allocateData(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo, DATA_GBODE *gbData)
 {
   DATA_GBODEF *gbfData = (DATA_GBODEF *)calloc(1, sizeof(DATA_GBODEF));
   gbData->gbfData = gbfData;
@@ -254,6 +254,7 @@ int gbodef_allocateData(DATA *data, threadData_t *threadData, DATA_GBODE *gbData
     infoStreamPrint(LOG_SOLVER, 0, "Linear interpolation is used for emitting results");
     break;
   case GB_INTERPOL_HERMITE:
+  case GB_INTERPOL_HERMITE_a:
   case GB_INTERPOL_HERMITE_b:
   case GB_INTERPOL_HERMITE_ERRCTRL:
     infoStreamPrint(LOG_SOLVER, 0, "Hermite interpolation is used for the slow states");
@@ -269,8 +270,11 @@ int gbodef_allocateData(DATA *data, threadData_t *threadData, DATA_GBODE *gbData
   if (ACTIVE_STREAM(LOG_GBODE_STATES))
   {
     char filename[4096];
-    sprintf(filename, "%s_ActiveStates.txt", data->modelData->modelFilePrefix);
+    unsigned int bufSize = 4096;
+    snprintf(filename, bufSize, "%s_ActiveStates.txt", data->modelData->modelFilePrefix);
     gbfData->fastStatesDebugFile = omc_fopen(filename, "w");
+    warningStreamPrint(LOG_STDOUT, 0, "LOG_GBODE_STATES sets -noEquidistantTimeGrid for emitting results!");
+    solverInfo->integratorSteps = TRUE;
   }
   else
   {
@@ -447,10 +451,11 @@ int gbode_allocateData(DATA *data, threadData_t *threadData, SOLVER_INFO *solver
   }
 
   char buffer[1024];
+  unsigned int bufSize = 1024;
   if (gbData->multi_rate) {
-    sprintf(buffer, "%s", " and slow states interpolation");
+    snprintf(buffer, bufSize, "%s", " and slow states interpolation");
   } else {
-    sprintf(buffer, "%s"," ");
+    snprintf(buffer, bufSize, "%s"," ");
   }
   switch (gbData->interpolation)
   {
@@ -458,6 +463,7 @@ int gbode_allocateData(DATA *data, threadData_t *threadData, SOLVER_INFO *solver
     infoStreamPrint(LOG_SOLVER, 0, "Linear interpolation is used for emitting results%s", buffer);
     break;
   case GB_INTERPOL_HERMITE_ERRCTRL:
+  case GB_INTERPOL_HERMITE_a:
   case GB_INTERPOL_HERMITE_b:
   case GB_INTERPOL_HERMITE:
     infoStreamPrint(LOG_SOLVER, 0, "Hermite interpolation is used for emitting results%s", buffer);
@@ -475,7 +481,7 @@ int gbode_allocateData(DATA *data, threadData_t *threadData, SOLVER_INFO *solver
 
 
   if (gbData->multi_rate) {
-    gbodef_allocateData(data, threadData, gbData);
+    gbodef_allocateData(data, threadData, solverInfo, gbData);
     gbData->tableau->isKRightAvailable = FALSE;
   } else {
     gbData->gbfData = NULL;
@@ -621,7 +627,7 @@ void gbodef_init(DATA* data, threadData_t* threadData, SOLVER_INFO* solverInfo)
   gbfData->didEventStep = FALSE;
 
   gbfData->time = gbData->time;
-  gbfData->stepSize = gbData->lastStepSize/2.5;
+  gbfData->stepSize = 0.1*gbData->stepSize*IController(&(gbData->err_fast), &(gbData->stepSize), 1);
 
   memcpy(gbfData->yOld, gbData->yOld, sizeof(double) * nStates);
   memcpy(gbfData->y, gbData->y, sizeof(double) * nStates);
@@ -833,11 +839,12 @@ int gbodef_main(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo, d
 
       // error handling: try half of the step size!
       if (integrator_step_info != 0) {
+        (gbfData->stats).nConvergenveTestFailures++;
         infoStreamPrint(LOG_SOLVER, 0, "gbodef_main: Failed to calculate step at time = %5g.", gbfData->time);
         gbfData->stepSize *= 0.5;
         infoStreamPrint(LOG_SOLVER, 0, "Try half of the step size = %g", gbfData->stepSize);
-        if (gbfData->stepSize < MINIMAL_STEP_SIZE) {
-          errorStreamPrint(LOG_STDOUT, 0, "Simulation aborted! Minimum step size %g reached, but error still to large.", MINIMAL_STEP_SIZE);
+        if (gbfData->stepSize < GB_MINIMAL_STEP_SIZE) {
+          errorStreamPrint(LOG_STDOUT, 0, "Simulation aborted! Minimum step size %g reached, but error still to large.", GB_MINIMAL_STEP_SIZE);
           messageClose(LOG_SOLVER);
           return -1;
         }
@@ -888,7 +895,7 @@ int gbodef_main(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo, d
         infoStreamPrint(LOG_SOLVER, 0, "Reject step from %10g to %10g, error %10g, new stepsize %10g",
                         gbfData->time, gbfData->time + gbfData->lastStepSize, err, gbfData->stepSize);
         if (ACTIVE_STREAM(LOG_GBODE_STATES)) {
-          dumpFastStates_gbf(gbData, gbfData->time + gbfData->lastStepSize);
+          dumpFastStates_gbf(gbData, gbfData->time + gbfData->lastStepSize, 1);
         }
       }
     } while (err > 1);
@@ -949,7 +956,7 @@ int gbodef_main(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo, d
       }
 
       if (ACTIVE_STREAM(LOG_GBODE_STATES)) {
-        dumpFastStates_gb(gbData, TRUE, eventTime);
+        dumpFastStates_gb(gbData, TRUE, eventTime, 0);
       }
 
       // Get out of the integration routine for event handling
@@ -986,7 +993,7 @@ int gbodef_main(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo, d
                     gbfData->time - gbfData->lastStepSize, gbfData->time, err, gbfData->stepSize);
 
     if (ACTIVE_STREAM(LOG_GBODE_STATES)) {
-      dumpFastStates_gbf(gbData, gbfData->time);
+      dumpFastStates_gbf(gbData, gbfData->time, 0);
     }
 
     /* emit step, if integratorSteps is selected */
@@ -1008,7 +1015,7 @@ int gbodef_main(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo, d
       }
     }
 
-    if ((gbData->timeRight - gbfData->time) < MINIMAL_STEP_SIZE || gbData->stepSize < MINIMAL_STEP_SIZE) {
+    if ((gbData->timeRight - gbfData->time) < GB_MINIMAL_STEP_SIZE || gbData->stepSize < GB_MINIMAL_STEP_SIZE) {
       gbfData->time = gbData->timeRight;
       break;
     }
@@ -1111,7 +1118,7 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
     // run multirate step
     gb_step_info = gbodef_main(data, threadData, solverInfo, targetTime);
     // synchronize y, yRight , kRight and buffer
-    if (fabs(gbData->timeRight - gbData->gbfData->timeRight) < MINIMAL_STEP_SIZE) {
+    if (fabs(gbData->timeRight - gbData->gbfData->timeRight) < GB_MINIMAL_STEP_SIZE) {
       gbData->time = gbData->timeRight;
       memcpy(gbData->y, gbData->gbfData->y, nStates * sizeof(double));
       memcpy(gbData->yOld, gbData->y, nStates * sizeof(double));
@@ -1130,7 +1137,7 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
 
       if (ACTIVE_STREAM(LOG_GBODE_STATES)) {
         // dump fast states in file
-        dumpFastStates_gb(gbData, FALSE, gbData->time);
+        dumpFastStates_gb(gbData, FALSE, gbData->time, 0);
       }
     }
     if (gb_step_info !=0) {
@@ -1188,20 +1195,29 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
 
       // error handling: try half of the step size!
       if (gb_step_info != 0) {
+        gbData->stats.nConvergenveTestFailures++;
         infoStreamPrint(LOG_SOLVER, 0, "gbode_main: Failed to calculate step at time = %5g.", gbData->time + gbData->stepSize);
         if (gbData->ctrl_method == GB_CTRL_CNST) {
           errorStreamPrint(LOG_STDOUT, 0, "Simulation aborted since gbode is running with fixed step size!");
           messageClose(LOG_SOLVER);
           return -1;
         } else {
-          if (gbData->stepSize > MINIMAL_STEP_SIZE) {
+          if (ACTIVE_STREAM(LOG_GBODE_STATES)) {
+            gbData->err_slow = 0;
+            gbData->err_fast = 0;
+            gbData->err_int = 0;
+            // dump fast states in file
+            dumpFastStates_gb(gbData, FALSE, gbData->time + gbData->stepSize, 3);
+          }
+
+          if (gbData->stepSize > GB_MINIMAL_STEP_SIZE) {
             // Try smaller steps, if possible.
             gbData->stepSize = gbData->stepSize / 2.;
             warningStreamPrint(LOG_SOLVER, 0, "Try half of the step size = %g", gbData->stepSize);
             err = 100;
             continue;
           } else {
-            errorStreamPrint(LOG_STDOUT, 0, "Simulation aborted because the step size is less then %g!", MINIMAL_STEP_SIZE);
+            errorStreamPrint(LOG_STDOUT, 0, "Simulation aborted because the step size is less then %g!", GB_MINIMAL_STEP_SIZE);
             messageClose(LOG_SOLVER);
             return -1;
           }
@@ -1285,7 +1301,7 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
         if (ACTIVE_STREAM(LOG_GBODE_STATES)) {
           // dump fast states in file
           gbData->err_slow = err;
-          dumpFastStates_gb(gbData, FALSE, gbData->time + gbData->lastStepSize);
+          dumpFastStates_gb(gbData, FALSE, gbData->time + gbData->lastStepSize, 1);
         }
         continue;
       }
@@ -1313,17 +1329,25 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
         printVector_gbf(LOG_GBODE_V, "e", gbData->errest, nStates, (gbData->timeLeft + gbData->timeRight)/2, gbData->nSlowStates, gbData->slowStatesIdx);
         messageClose(LOG_GBODE_V);
       }
+      if (gbData->ctrl_method != GB_CTRL_CNST && ((gbData->interpolation == GB_INTERPOL_HERMITE_ERRCTRL)  || (gbData->interpolation == GB_DENSE_OUTPUT_ERRCTRL))) {
+        gbData->errValues[0] = fmax(err, gbData->err_int);
+        gbData->stepSize = gbData->lastStepSize * gbData->stepSize_control(gbData->errValues, gbData->stepSizeValues, gbData->tableau->error_order);
+      }
       // reject step, if interpolaton error is too large
       if (( gbData->nFastStates>0) && (gbData->err_int > 1 ) && gbData->ctrl_method != GB_CTRL_CNST &&
           ((gbData->interpolation == GB_INTERPOL_HERMITE_ERRCTRL)  || (gbData->interpolation == GB_DENSE_OUTPUT_ERRCTRL))) {
         err = 100;
-        gbData->stepSize = gbData->lastStepSize*IController(&(gbData->err_int), &(gbData->lastStepSize), 1);
+        if (gbData->stepSize < GB_MINIMAL_STEP_SIZE) {
+          errorStreamPrint(LOG_STDOUT, 0, "Simulation aborted! Minimum step size %g reached, but interpolation error still to large.", GB_MINIMAL_STEP_SIZE);
+          messageClose(LOG_SOLVER);
+          return -1;
+        }
         infoStreamPrint(LOG_SOLVER, 0, "Reject step from %10g to %10g, interpolation error %10g, new stepsize %10g",
                         gbData->time, gbData->time + gbData->lastStepSize, gbData->err_int, gbData->stepSize);
 
         if (ACTIVE_STREAM(LOG_GBODE_STATES)) {
           // dump fast states in file
-          dumpFastStates_gb(gbData, FALSE, gbData->time + gbData->lastStepSize);
+          dumpFastStates_gb(gbData, FALSE, gbData->time + gbData->lastStepSize, 2);
         }
 
         // count failed steps and output information on the solver status
@@ -1347,10 +1371,14 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
           printIntVector_gb(LOG_GBODE, "sr", gbData->sortedStatesIdx, nStates, gbData->timeRight);
           messageClose(LOG_GBODE);
         }
+        if (ACTIVE_STREAM(LOG_GBODE_STATES)) {
+          // dump fast states in file
+          dumpFastStates_gb(gbData, FALSE, gbData->time + gbData->lastStepSize, -1);
+        }
         // run multirate step
         gb_step_info = gbodef_main(data, threadData, solverInfo, targetTime);
         // synchronize relevant information
-        if (fabs(gbData->timeRight - gbData->gbfData->timeRight) < MINIMAL_STEP_SIZE) {
+        if (fabs(gbData->timeRight - gbData->gbfData->timeRight) < GB_MINIMAL_STEP_SIZE) {
           memcpy(gbData->y, gbData->gbfData->y, nStates * sizeof(double));
           memcpy(gbData->yRight, gbData->gbfData->yRight, nStates * sizeof(double));
           memcpy(gbData->kRight, gbData->gbfData->kRight, nStates * sizeof(double));
@@ -1415,7 +1443,7 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
 
         if (ACTIVE_STREAM(LOG_GBODE_STATES)) {
           // dump fast states in file
-          dumpFastStates_gb(gbData, TRUE, eventTime);
+          dumpFastStates_gb(gbData, TRUE, eventTime, 0);
         }
 
         // return to solver main routine for proper event handling (iteration)
@@ -1445,7 +1473,7 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
 
     if (ACTIVE_STREAM(LOG_GBODE_STATES)) {
       // dump fast states in file
-      dumpFastStates_gb(gbData, FALSE, gbData->time);
+      dumpFastStates_gb(gbData, FALSE, gbData->time, 0);
     }
 
     /* emit step, if integratorSteps is selected */
@@ -1467,7 +1495,7 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
       }
     }
 
-    if ((stopTime - gbData->time) < MINIMAL_STEP_SIZE)
+    if ((stopTime - gbData->time) < GB_MINIMAL_STEP_SIZE)
     {
       gbData->time = stopTime;
       break;
@@ -1517,10 +1545,11 @@ int gbode_birate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverInfo)
   if (!gbData->isExplicit)
     gbData->stats.nCallsJacobian = gbData->nlsData->numberOfJEval;
 
-  if (fabs(targetTime - stopTime) < MINIMAL_STEP_SIZE && ACTIVE_STREAM(LOG_STATS)) {
+  if (fabs(targetTime - stopTime) < GB_MINIMAL_STEP_SIZE && ACTIVE_STREAM(LOG_STATS)) {
     infoStreamPrint(LOG_STATS, 0, "gbode (birate integration): slow: %s / fast: %s",
                     GB_METHOD_NAME[gbData->GM_method], GB_METHOD_NAME[gbData->gbfData->GM_method]);
     logSolverStats(LOG_STATS, "inner integration", stopTime, stopTime, 0, &gbData->gbfData->stats);
+    logSolverStats(LOG_STATS, "outer integration", stopTime, stopTime, 0, &gbData->stats);
   }
   /* Write statistics to the solverInfo data structure */
   logSolverStats(LOG_SOLVER_V, "gb_singlerate", solverInfo->currentTime, gbData->time, gbData->stepSize, &gbData->stats);
@@ -1641,6 +1670,7 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn
 
       // error handling: try half of the step size!
       if (gb_step_info != 0) {
+        gbData->stats.nConvergenveTestFailures++;
         infoStreamPrint(LOG_SOLVER, 0, "gbode_main: Failed to calculate step at time = %5g.", gbData->time + gbData->stepSize);
         if (gbData->ctrl_method == GB_CTRL_CNST) {
           errorStreamPrint(LOG_STDOUT, 0, "Simulation aborted since gbode is running with fixed step size!");
@@ -1649,8 +1679,8 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn
         } else {
           gbData->stepSize *= 0.5;
           infoStreamPrint(LOG_SOLVER, 0, "Try half of the step size = %g", gbData->stepSize);
-          if (gbData->stepSize < MINIMAL_STEP_SIZE) {
-            errorStreamPrint(LOG_STDOUT, 0, "Simulation aborted! Minimum step size %g reached, but error still to large.", MINIMAL_STEP_SIZE);
+          if (gbData->stepSize < GB_MINIMAL_STEP_SIZE) {
+            errorStreamPrint(LOG_STDOUT, 0, "Simulation aborted! Minimum step size %g reached, but error still to large.", GB_MINIMAL_STEP_SIZE);
             messageClose(LOG_SOLVER);
             return -1;
           }
@@ -1700,7 +1730,7 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn
       }
       memcpy(gbData->kRight, fODE, nStates * sizeof(double));
 
-      if (gbData->ctrl_method != GB_CTRL_CNST && ((gbData->interpolation == GB_INTERPOL_HERMITE_ERRCTRL)  || (gbData->interpolation == GB_DENSE_OUTPUT_ERRCTRL))) {
+      if (ACTIVE_STREAM(LOG_SOLVER) || (gbData->ctrl_method != GB_CTRL_CNST && ((gbData->interpolation == GB_INTERPOL_HERMITE_ERRCTRL)  || (gbData->interpolation == GB_DENSE_OUTPUT_ERRCTRL)))) {
         gbData->err_int = error_interpolation_gb(gbData, nStates, NULL, Rtol);
       }
       if (ACTIVE_STREAM(LOG_GBODE_V)) {
@@ -1713,11 +1743,20 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn
         printVector_gbf(LOG_GBODE_V, "e", gbData->errest, nStates, (gbData->timeLeft + gbData->timeRight)/2, gbData->nSlowStates, gbData->slowStatesIdx);
         messageClose(LOG_GBODE_V);
       }
+      if (gbData->ctrl_method != GB_CTRL_CNST && ((gbData->interpolation == GB_INTERPOL_HERMITE_ERRCTRL)  || (gbData->interpolation == GB_DENSE_OUTPUT_ERRCTRL))) {
+        gbData->errValues[0] = fmax(err, gbData->err_int);
+        gbData->stepSize = gbData->lastStepSize * gbData->stepSize_control(gbData->errValues, gbData->stepSizeValues, gbData->tableau->error_order);
+      }
       // reject step, if interpolaton error is too large
       if ((gbData->err_int > 1 ) && gbData->ctrl_method != GB_CTRL_CNST &&
           ((gbData->interpolation == GB_INTERPOL_HERMITE_ERRCTRL)  || (gbData->interpolation == GB_DENSE_OUTPUT_ERRCTRL))) {
         err = 100;
-        gbData->stepSize = gbData->lastStepSize*IController(&(gbData->err_int), &(gbData->lastStepSize), 1);
+        // gbData->stepSize = gbData->lastStepSize*IController(&(gbData->err_int), &(gbData->lastStepSize), 1);
+        if (gbData->stepSize < GB_MINIMAL_STEP_SIZE) {
+          errorStreamPrint(LOG_STDOUT, 0, "Simulation aborted! Minimum step size %g reached, but interpolation error still to large.", GB_MINIMAL_STEP_SIZE);
+          messageClose(LOG_SOLVER);
+          return -1;
+        }
         infoStreamPrint(LOG_SOLVER, 0, "Reject step from %10g to %10g, interpolation error %10g, new stepsize %10g",
                         gbData->time, gbData->time + gbData->lastStepSize, gbData->err_int, gbData->stepSize);
 
@@ -1752,25 +1791,48 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn
     // check for events, if event is detected stop integrator and trigger event iteration
     eventTime = checkForEvents(data, threadData, solverInfo, gbData->timeLeft, gbData->yLeft, gbData->timeRight, gbData->yRight, FALSE, &foundEvent);
     if (foundEvent) {
-      solverInfo->currentTime = eventTime;
-      sData->timeValue = eventTime;
+      if (eventTime < targetTime + solverInfo->currentStepSize/2)
+      {
+        solverInfo->currentTime = eventTime;
+        sData->timeValue = eventTime;
 
-      // sData->realVars are the "numerical" values on the right hand side of the event (hopefully)
-      gbData->time = eventTime;
-      memcpy(gbData->yOld, sData->realVars, gbData->nStates * sizeof(double));
+        // sData->realVars are the "numerical" values on the right hand side of the event (hopefully)
+        gbData->time = eventTime;
+        memcpy(gbData->yOld, sData->realVars, gbData->nStates * sizeof(double));
 
-      /* write statistics to the solverInfo data structure */
-      setSolverStats(solverInfo->solverStatsTmp, &gbData->stats);
+        /* write statistics to the solverInfo data structure */
+        setSolverStats(solverInfo->solverStatsTmp, &gbData->stats);
+
+        // log the emitted result
+        if (ACTIVE_STREAM(LOG_GBODE)){
+          infoStreamPrint(LOG_GBODE, 1, "Emit result (single-rate integration):");
+          printVector_gb(LOG_GBODE, " y", sData->realVars, nStates, sData->timeValue);
+          messageClose(LOG_GBODE);
+        }
+        // return to solver main routine for proper event handling (iteration)
+        messageClose(LOG_SOLVER);
+        return 0;
+      } else {
+        // ToDo: If the solver does large steps and finds an event, the interpolation is
+        // done in solver_main (linearly) and therefore the states are not very well approximated.
+        // Current solution: Step back to the communication interval before the event and event detection
+        // needs to be repeated
+        listClear(solverInfo->eventLst);
+        gbData->lastStepSize = (eventTime - solverInfo->currentStepSize/2) - gbData->time;
+        sData->timeValue = (eventTime - solverInfo->currentStepSize/2);
+        gb_interpolation(gbData->interpolation,
+                        gbData->timeLeft,  gbData->yLeft,  gbData->kLeft,
+                        gbData->timeRight, gbData->yRight, gbData->kRight,
+                                sData->timeValue,  sData->realVars,
+                        nStates, NULL, nStates, gbData->tableau, gbData->x, gbData->k);
+        memcpy(gbData->y, sData->realVars, gbData->nStates * sizeof(double));
+
+        gbData->timeRight = sData->timeValue;
+        memcpy(gbData->yRight, sData->realVars, gbData->nStates * sizeof(double));
+        gbode_fODE(data, threadData, &(gbData->stats.nCallsODE));
+        memcpy(gbData->kRight, fODE, nStates * sizeof(double));
 
-      // log the emitted result
-      if (ACTIVE_STREAM(LOG_GBODE)){
-        infoStreamPrint(LOG_GBODE, 1, "Emit result (single-rate integration):");
-        printVector_gb(LOG_GBODE, " y", sData->realVars, nStates, sData->timeValue);
-        messageClose(LOG_GBODE);
       }
-      // return to solver main routine for proper event handling (iteration)
-      messageClose(LOG_SOLVER);
-      return 0;
     }
 
     /* update time with performed stepSize */
@@ -1779,7 +1841,7 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn
     /* step is accepted and yOld needs to be updated */
     memcpy(gbData->yOld, gbData->y, nStates * sizeof(double));
     infoStreamPrint(LOG_SOLVER, 0, "Accept step from %10g to %10g, error %10g interpolation error %10g, new stepsize %10g",
-                    gbData->time - gbData->lastStepSize, gbData->time, gbData->errValues[0], gbData->err_int, gbData->stepSize);
+                    gbData->time - gbData->lastStepSize, gbData->time, err, gbData->err_int, gbData->stepSize);
 
     // Rotate buffer
     for (i = (gbData->ringBufferSize - 1); i > 0 ; i--) {
@@ -1812,7 +1874,7 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn
     }
 
     // stop, if simulation nearly reached stopTime
-    if ((stopTime - gbData->time) < MINIMAL_STEP_SIZE) {
+    if ((stopTime - gbData->time) < GB_MINIMAL_STEP_SIZE) {
       gbData->time = stopTime;
       break;
     }
@@ -1853,7 +1915,7 @@ int gbode_singlerate(DATA *data, threadData_t *threadData, SOLVER_INFO *solverIn
   /* Solver statistics */
   if (!gbData->isExplicit)
     gbData->stats.nCallsJacobian = gbData->nlsData->numberOfJEval;
-  if (fabs(targetTime - stopTime) < MINIMAL_STEP_SIZE && ACTIVE_STREAM(LOG_STATS)) {
+  if (fabs(targetTime - stopTime) < GB_MINIMAL_STEP_SIZE && ACTIVE_STREAM(LOG_STATS)) {
     infoStreamPrint(LOG_STATS, 0, "gbode (single-rate integration): %s", GB_METHOD_NAME[gbData->GM_method]);
   }
   /* Write statistics to the solverInfo data structure */
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_nls.c b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_nls.c
index 845ae6bf03b..2d2506191bf 100644
--- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_nls.c
+++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_nls.c
@@ -46,7 +46,21 @@
 #include "nonlinearSystem.h"
 
 #include "util/jacobian_util.h"
+#include "util/rtclock.h"
 
+/**
+ * @brief Specific error handling of kinsol for gbode
+ *
+ * @param error_code  Reported error code
+ * @param module      Module of failure
+ * @param function    Nonlinear function
+ * @param msg         Message of failure
+ * @param data        Pointer to userData
+ */
+void GB_KINErrHandler(int error_code, const char *module, const char *function, char *msg, void *data) {
+// Preparation for specific error handling of the solution process of kinsol for gbode
+// This is needed to speed up simulation in case of failure
+}
 
 /**
  * @brief Initialize static data of non-linear system for DIRK.
@@ -264,8 +278,14 @@ NONLINEAR_SYSTEM_DATA* initRK_NLS_DATA(DATA* data, threadData_t* threadData, DAT
     solverData->initHomotopyData = NULL;
     nlsData->solverData = solverData;
 
-    int flag = KINSetNumMaxIters(((NLS_KINSOL_DATA*)solverData->ordinaryData)->kinsolMemory, nlsData->size * 10);
+    int flag;
+    NLS_KINSOL_DATA* kin_mem = ((NLS_KINSOL_DATA*)solverData->ordinaryData)->kinsolMemory;
+    flag = KINSetNumMaxIters(kin_mem, nlsData->size * 4);
     checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNumMaxIters");
+    flag = KINSetMaxSetupCalls(kin_mem, 10);
+    checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetMaxSetupCalls");
+    flag = KINSetErrHandlerFn(kin_mem, GB_KINErrHandler, NULL);
+    checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetErrHandlerFn");
     break;
   default:
     throwStreamPrint(NULL, "Memory allocation for NLS method %s not yet implemented.", GB_NLS_METHOD_NAME[gbData->nlsSolverMethod]);
@@ -389,6 +409,96 @@ void freeRK_NLS_DATA(NONLINEAR_SYSTEM_DATA* nlsData) {
   return;
 }
 
+/**
+ * @brief Set kinsol parameters
+ *
+ * @param kin_mem       Pointer to kinsol data object
+ * @param numIter       Number of nonlinear iterations
+ * @param jacUpdate     Update of jacobian necessary (SUNFALSE => yes)
+ * @param maxJacUpdate  Maximal number of constant jacobian
+ */
+void set_kinsol_parameters(NLS_KINSOL_DATA* kin_mem, int numIter, int jacUpdate, int maxJacUpdate) {
+    int flag;
+
+    flag = KINSetNumMaxIters(kin_mem, numIter);
+    checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNumMaxIters");
+    flag = KINSetNoInitSetup(kin_mem, jacUpdate);
+    checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetNoInitSetup");
+    flag = KINSetMaxSetupCalls(kin_mem, maxJacUpdate);
+    checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINSetMaxSetupCalls");
+
+}
+
+/**
+ * @brief Get the kinsol statistics object
+ *
+ * @param kin_mem Pointer to kinsol data object
+ */
+void get_kinsol_statistics(NLS_KINSOL_DATA* kin_mem) {
+  int flag;
+  long int nIters, nFuncEvals, nJacEvals;
+  double fnorm;
+
+  // Get number of nonlinear iteration steps
+  flag = KINGetNumNonlinSolvIters(kin_mem, &nIters);
+  checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINGetNumNonlinSolvIters");
+
+  // Get the error of the residual function
+  flag = KINGetFuncNorm(kin_mem, &fnorm);
+  checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINGetFuncNorm");
+
+  // Get the number of jacobian evaluation
+  flag = KINGetNumJacEvals(kin_mem, &nJacEvals);
+  checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINGetNumJacEvals");
+
+  // Get the number of function evaluation
+  flag = KINGetNumFuncEvals(kin_mem, &nFuncEvals);
+  checkReturnFlag_SUNDIALS(flag, SUNDIALS_KIN_FLAG, "KINGetNumFuncEvals");
+
+  // Report numbers
+  infoStreamPrint(LOG_GBODE_NLS, 0, "Kinsol statistics: nIters = %ld, nFuncEvals = %ld, nJacEvals = %ld,  fnorm:  %14.12g", nIters, nFuncEvals, nJacEvals, fnorm);
+}
+/**
+ * @brief Special treatment when solving non linear systems of equations
+ *
+ *        Will be described, when it is ready
+ *
+ * @param data                Pointer to runtime data struct.
+ * @param threadData          Thread data for error handling.
+ * @param nlsData             Non-linear solver data.
+ * @param gbData              Runge-Kutta method.
+ * @return NLS_SOLVER_STATUS  Return NLS_SOLVED on success and NLS_FAILED otherwise.
+ */
+NLS_SOLVER_STATUS solveNLS_gb(DATA *data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA* nlsData, DATA_GBODE* gbData) {
+  struct dataSolver * solverData = (struct dataSolver *)nlsData->solverData;
+  NLS_SOLVER_STATUS solved;
+
+  // Debug nonlinear solution process
+  rtclock_t clock;
+  double cpu_time_used;
+
+  if (ACTIVE_STREAM(LOG_GBODE_NLS)) {
+    rt_ext_tp_tick(&clock);
+  }
+
+  if (gbData->nlsSolverMethod == GB_NLS_KINSOL) {
+    // Get kinsol data object
+    NLS_KINSOL_DATA* kin_mem = ((NLS_KINSOL_DATA*)solverData->ordinaryData)->kinsolMemory;
+
+    set_kinsol_parameters(kin_mem, nlsData->size * 4, SUNFALSE, 10);
+    solved = solveNLS(data, threadData, nlsData);
+    if (ACTIVE_STREAM(LOG_GBODE_NLS)) get_kinsol_statistics(kin_mem);
+  } else {
+    solved = solveNLS(data, threadData, nlsData);
+  }
+
+  if (ACTIVE_STREAM(LOG_GBODE_NLS)) {
+      cpu_time_used = rt_ext_tp_tock(&clock);
+      infoStreamPrint(LOG_GBODE_NLS, 0, "Time needed for solving the NLS:  %20.16g", cpu_time_used);
+  }
+
+  return solved;
+}
 
 /**
  * @brief Residual function for non-linear system of generic multi-step methods.
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_nls.h b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_nls.h
index 87662d61099..ae0135cd2ff 100644
--- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_nls.h
+++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_nls.h
@@ -49,6 +49,9 @@ NONLINEAR_SYSTEM_DATA* initRK_NLS_DATA(DATA* data, threadData_t* threadData, DAT
 NONLINEAR_SYSTEM_DATA* initRK_NLS_DATA_MR(DATA* data, threadData_t* threadData, DATA_GBODEF* gbfData);
 void freeRK_NLS_DATA( NONLINEAR_SYSTEM_DATA* nlsData);
 
+//Specific treatment of NLS within gbode
+NLS_SOLVER_STATUS solveNLS_gb(DATA *data, threadData_t *threadData, NONLINEAR_SYSTEM_DATA* nonlinsys, DATA_GBODE* gbData);
+
 // Residuum and Jacobian functions for diagonal implicit (DIRK) and implicit (IRK) Runge-Kutta methods.
 void residual_MS(RESIDUAL_USERDATA* userData, const double *xloc, double *res, const int *iflag);
 void residual_DIRK(RESIDUAL_USERDATA* userData, const double *xloc, double *res, const int *iflag);
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_step.c b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_step.c
index 0e16bb70ce4..33faf1caca2 100644
--- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_step.c
+++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_step.c
@@ -32,8 +32,8 @@
  */
 
 #include "gbode_main.h"
+#include "gbode_nls.h"
 #include "gbode_util.h"
-#include "util/rtclock.h"
 
 /**
  * @brief Generic multi-step function.
@@ -95,18 +95,7 @@ int full_implicit_MS(DATA* data, threadData_t* threadData, SOLVER_INFO* solverIn
   memcpy(nlsData->nlsxOld, nlsData->nlsx, nStates*sizeof(modelica_real));
   memcpy(nlsData->nlsxExtrapolation, nlsData->nlsx, nStates*sizeof(modelica_real));
 
-  if (ACTIVE_STREAM(LOG_GBODE_NLS)) {
-    rtclock_t clock;
-    double cpu_time_used;
-
-    rt_ext_tp_tick(&clock);
-    solved = solveNLS(data, threadData, nlsData);
-    cpu_time_used = rt_ext_tp_tock(&clock);
-
-    infoStreamPrint(LOG_GBODE_NLS, 0, "time needed for a solving NLS:  %20.16g", cpu_time_used);
-  } else {
-    solved = solveNLS(data, threadData, nlsData);
-  }
+  solved = solveNLS_gb(data, threadData, nlsData, gbData);
 
   if (solved != NLS_SOLVED) {
     warningStreamPrint(LOG_SOLVER, 0, "gbode error: Failed to solve NLS in full_implicit_MS");
@@ -201,18 +190,7 @@ int full_implicit_MS_MR(DATA* data, threadData_t* threadData, SOLVER_INFO* solve
   memcpy(nlsData->nlsxOld, nlsData->nlsx, nStates*sizeof(modelica_real));
   memcpy(nlsData->nlsxExtrapolation, nlsData->nlsx, nStates*sizeof(modelica_real));
 
-  if (ACTIVE_STREAM(LOG_GBODE_NLS_V)) {
-    rtclock_t clock;
-    double cpu_time_used;
-
-    rt_ext_tp_tick(&clock);
-    solved = solveNLS(data, threadData, nlsData);
-    cpu_time_used = rt_ext_tp_tock(&clock);
-
-    infoStreamPrint(LOG_GBODE_NLS_V, 0, "time needed for a solving NLS:  %20.16g", cpu_time_used);
-  } else {
-    solved = solveNLS(data, threadData, nlsData);
-  }
+  solved = solveNLS_gb(data, threadData, nlsData, gbData);
 
   if (solved != NLS_SOLVED) {
     warningStreamPrint(LOG_SOLVER, 0, "gbodef error: Failed to solve NLS in full_implicit_MS_MR");
@@ -309,40 +287,29 @@ int expl_diag_impl_RK(DATA* data, threadData_t* threadData, SOLVER_INFO* solverI
     } else {
       // solve for x: 0 = yold-x + h*(sum(A[i,j]*k[j], i=1..j-1) + A[i,i]*f(t + c[i]*h, x))
       NONLINEAR_SYSTEM_DATA* nlsData = gbData->nlsData;
+      struct dataSolver * solverData = (struct dataSolver *)nlsData->solverData;
+      NLS_KINSOL_DATA* kin_mem = ((NLS_KINSOL_DATA*)solverData->ordinaryData)->kinsolMemory;
 
       // Set start vector
       memcpy(nlsData->nlsx,    gbData->yOld, nStates*sizeof(modelica_real));
       memcpy(nlsData->nlsxOld, gbData->yOld, nStates*sizeof(modelica_real));
       extrapolation_gb(gbData, nlsData->nlsxExtrapolation, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
 
-      // Debug nonlinear solution process
-      if (ACTIVE_STREAM(LOG_GBODE_NLS)) {
-        rtclock_t clock;
-        double cpu_time_used;
-
-        //Solve nonlinear equation system
-        rt_ext_tp_tick(&clock);
-        solved = solveNLS(data, threadData, nlsData);
-        cpu_time_used = rt_ext_tp_tock(&clock);
+      solved = solveNLS_gb(data, threadData, nlsData, gbData);
 
-        infoStreamPrint(LOG_GBODE_NLS, 0, "time needed for solving the NLS:  %20.16g", cpu_time_used);
-      } else {
-        //Solve nonlinear equation system
-        solved = solveNLS(data, threadData, nlsData);
+      if (solved != NLS_SOLVED) {
+        warningStreamPrint(LOG_SOLVER, 0, "gbode error: Failed to solve NLS in expl_diag_impl_RK in stage %d", stage_);
+        return -1;
       }
 
       if (ACTIVE_STREAM(LOG_GBODE_NLS_V)) {
         infoStreamPrint(LOG_GBODE_NLS_V, 1, "NLS - start values and solution of the NLS:");
-        printVector_gb(LOG_GBODE_NLS_V, "xo", nlsData->nlsxOld, nStates, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
+        printVector_gb(LOG_GBODE_NLS_V, "x0", nlsData->nlsxOld, nStates, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
         printVector_gb(LOG_GBODE_NLS_V, "xS", nlsData->nlsxExtrapolation, nStates, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
         printVector_gb(LOG_GBODE_NLS_V, "xL", nlsData->nlsx,              nStates, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
         messageClose(LOG_GBODE_NLS_V);
       }
 
-      if (solved != NLS_SOLVED) {
-        warningStreamPrint(LOG_SOLVER, 0, "gbode error: Failed to solve NLS in expl_diag_impl_RK in stage %d", stage_);
-        return -1;
-      }
       memcpy(gbData->x + stage_ * nStates, nlsData->nlsx, nStates*sizeof(double));
     }
     // copy last calculation of fODE, which should coincide with k[i], here, it yields stage == stage_
@@ -457,18 +424,11 @@ int expl_diag_impl_RK_MR(DATA* data, threadData_t* threadData, SOLVER_INFO* solv
       extrapolation_gbf(gbData, gbData->y1, gbfData->time + gbfData->tableau->c[stage_] * gbfData->stepSize);
       projVector_gbf(nlsData->nlsxExtrapolation, gbData->y1, nFastStates, gbData->fastStatesIdx);
 
-      // Solve corresponding NLS
-      if (ACTIVE_STREAM(LOG_GBODE_NLS_V)) {
-        rtclock_t clock;
-        double cpu_time_used;
-
-        rt_ext_tp_tick(&clock);
-        solved = solveNLS(data, threadData, nlsData);
-        cpu_time_used = rt_ext_tp_tock(&clock);
+      solved = solveNLS_gb(data, threadData, nlsData, gbData);
 
-        infoStreamPrint(LOG_GBODE_NLS_V, 0, "time needed for a solving NLS:  %20.16g", cpu_time_used);
-      } else {
-        solved = solveNLS(data, threadData, nlsData);
+      if (solved != NLS_SOLVED) {
+        warningStreamPrint(LOG_SOLVER, 0, "gbodef error: Failed to solve NLS in expl_diag_impl_RK_MR in stage %d", stage_);
+        return -1;
       }
 
       // debug residuals
@@ -478,11 +438,6 @@ int expl_diag_impl_RK_MR(DATA* data, threadData_t* threadData, SOLVER_INFO* solv
         printVector_gb(LOG_GBODE_NLS, "xL", nlsData->nlsx,              nFastStates, gbfData->time + gbfData->tableau->c[stage_] * gbfData->stepSize);
         messageClose(LOG_GBODE_NLS);
       }
-
-      if (solved != NLS_SOLVED) {
-        warningStreamPrint(LOG_SOLVER, 0, "gbodef error: Failed to solve NLS in expl_diag_impl_RK_MR in stage %d", stage_);
-        return -1;
-      }
     }
 
     // copy last values of sData->realVars and fODE, which should coincide with x[i] and k[i]
@@ -552,17 +507,12 @@ int full_implicit_RK(DATA* data, threadData_t* threadData, SOLVER_INFO* solverIn
     extrapolation_gb(gbData, nlsData->nlsxExtrapolation + stage_*nStates, gbData->time + gbData->tableau->c[stage_] * gbData->stepSize);
   }
 
-  if (ACTIVE_STREAM(LOG_GBODE_NLS)) {
-    rtclock_t clock;
-    double cpu_time_used;
-
-    rt_ext_tp_tick(&clock);
-    solved = solveNLS(data, threadData, nlsData);
-    cpu_time_used = rt_ext_tp_tock(&clock);
+  solved = solveNLS_gb(data, threadData, nlsData, gbData);
 
-    infoStreamPrint(LOG_GBODE_NLS, 0, "time needed for a solving NLS:  %20.16g", cpu_time_used);
-  } else {
-    solved = solveNLS(data, threadData, nlsData);
+  if (solved != NLS_SOLVED) {
+    gbData->stats.nConvergenveTestFailures++;
+    warningStreamPrint(LOG_SOLVER, 0, "gbode error: Failed to solve NLS in full_implicit_RK");
+    return -1;
   }
 
   if (ACTIVE_STREAM(LOG_GBODE_NLS)) {
@@ -574,11 +524,6 @@ int full_implicit_RK(DATA* data, threadData_t* threadData, SOLVER_INFO* solverIn
     messageClose(LOG_GBODE_NLS);
   }
 
-  if (solved != NLS_SOLVED) {
-    gbData->stats.nConvergenveTestFailures++;
-    warningStreamPrint(LOG_SOLVER, 0, "gbode error: Failed to solve NLS in full_implicit_RK");
-    return -1;
-  }
 
   // Apply RK-scheme for determining the approximations at (gbData->time + gbData->stepSize)
   // y       = yold+h*sum(b[stage_]  * k[stage_], stage_=1..nStages);
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c
index 36d932c5cee..fd3c17caa79 100644
--- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c
+++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_tableau.c
@@ -1238,7 +1238,7 @@ void printButcherTableau(BUTCHER_TABLEAU* tableau) {
     }
     ct = snprintf(buffer, buffSize, "%s | ", line);
       for (j = 0; j<tableau->nStages; j++) {
-        ct += snprintf(buffer+ct, buffSize, "%s", line);
+        ct += snprintf(buffer+ct, buffSize-ct, "%s", line);
       }
     infoStreamPrint(LOG_SOLVER, 0, "%s", buffer);
     ct = snprintf(buffer, buffSize, "%10s | ", "");
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_util.c b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_util.c
index 675eb2b4e83..522ffdd5400 100644
--- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_util.c
+++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_util.c
@@ -239,6 +239,58 @@ void hermite_interpolation_b(double ta, double* fa, double tb, double* fb, doubl
   return;
 }
 
+/**
+ * @brief Hermite interpolation of specific vector components (only left derivative used)
+ *
+ * @param ta      Time value at the left hand side
+ * @param fa      Function values at the left hand side
+ * @param dfa     Derivative function values at the left hand side
+ * @param tb      Time value at the right hand side
+ * @param fb      Function values at the right hand side
+ * @param t       Time value at the interpolated time point
+ * @param f       Function values at the interpolated time point
+ * @param n       Size of vector f or size of index vector if non-NULL.
+ * @param idx     Index vector, can be NULL.
+ *                Specifies which parts of f should be interpolated.
+ */
+void hermite_interpolation_a(double ta, double* fa, double* dfa, double tb, double* fb, double t, double* f, int n, int* idx)
+{
+  double tat,tbt,tbta, h00, h01, h10;
+  int i, ii;
+
+  // omit division by zero
+  if (fabs(tb-ta) <= GBODE_EPSILON) {
+    if(idx != NULL) {
+      copyVector_gbf(f, fb, n, idx);
+    } else {
+      memcpy(f, fb, n*sizeof(double));
+    }
+    return;
+  }
+
+  tat  = (ta-t);
+  tbt  = (tb-t);
+  tbta = (tb-ta);
+  h01  = tat*tat/(tbta*tbta);
+  h00  = 1 - h01;
+  h10  = -tat*tbt/tbta;
+
+  if (idx == NULL) {
+    for (i=0; i<n; i++)
+    {
+      f[i] = h00*fa[i]+h01*fb[i]+h10*dfa[i];
+    }
+  } else {
+    for (ii=0; ii<n; ii++)
+    {
+      i = idx[ii];
+      f[i] = h00*fa[i]+h01*fb[i]+h10*dfa[i];
+    }
+  }
+
+  return;
+}
+
 /**
  * @brief Hermite interpolation of specific vector components
  *
@@ -270,6 +322,9 @@ void gb_interpolation(enum GB_INTERPOL_METHOD interpolMethod, double ta, double*
       tableau->dense_output(tableau, fa, x, k, (t - ta)/(tb-ta), (tb - ta), f, nIdx, idx, nStates);
       break;
     }
+  case GB_INTERPOL_HERMITE_a:
+    hermite_interpolation_a(ta, fa, dfa, tb, fb, t, f, nIdx, idx);
+    break;
   case GB_INTERPOL_HERMITE_b:
     hermite_interpolation_b(ta, fa, tb, fb, dfb, t, f, nIdx, idx);
     break;
@@ -281,6 +336,7 @@ void gb_interpolation(enum GB_INTERPOL_METHOD interpolMethod, double ta, double*
     throwStreamPrint(NULL, "Not handled case in gb_interpolation. Unknown interpolation method %i.", interpolMethod);
   }
 }
+
 /**
  * @brief  Difference between linear and hermite interpolation at intermediate points.
  *
@@ -290,18 +346,17 @@ double error_interpolation_gb(DATA_GBODE* gbData, int nIdx, int* idx, double tol
   int i, ii;
   double errint = 0.0, errtol;
 
- if (gbData->interpolation == GB_INTERPOL_HERMITE_ERRCTRL || gbData->interpolation == GB_INTERPOL_HERMITE || gbData->interpolation == GB_INTERPOL_HERMITE_b ) {
-    linear_interpolation(gbData->timeLeft,  gbData->yLeft,
-                        gbData->timeRight, gbData->yRight,
-                        (gbData->timeLeft + gbData->timeRight)/2, gbData->y1,
-                         nIdx, idx);
- } else {
-   gb_interpolation(gbData->interpolation, gbData->timeLeft,  gbData->yLeft,  gbData->kLeft,
-                     gbData->timeRight, gbData->yRight, gbData->kRight,
-                     (gbData->timeLeft + gbData->timeRight)/2, gbData->y1,
-                      nIdx, idx, gbData->nStates, gbData->tableau, gbData->x, gbData->k);
-
- }
+  if (gbData->interpolation == GB_DENSE_OUTPUT_ERRCTRL || gbData->interpolation == GB_DENSE_OUTPUT) {
+    gb_interpolation(gbData->interpolation, gbData->timeLeft,  gbData->yLeft,  gbData->kLeft,
+                      gbData->timeRight, gbData->yRight, gbData->kRight,
+                      (gbData->timeLeft + gbData->timeRight)/2, gbData->y1,
+                        nIdx, idx, gbData->nStates, gbData->tableau, gbData->x, gbData->k);
+  } else {
+    hermite_interpolation_a(gbData->timeLeft,  gbData->yLeft,  gbData->kLeft,
+                            gbData->timeRight, gbData->yRight,
+                            (gbData->timeLeft + gbData->timeRight)/2, gbData->y1,
+                            nIdx, idx);
+  }
   hermite_interpolation(gbData->timeLeft,  gbData->yLeft,  gbData->kLeft,
                         gbData->timeRight, gbData->yRight, gbData->kRight,
                         (gbData->timeLeft + gbData->timeRight)/2, gbData->errest,
@@ -436,7 +491,9 @@ void debugRingBuffer(enum LOG_STREAM stream, double* x, double* k, int nStates,
 }
 
 /**
- * @brief Prints a vector
+ * @brief Prints a vector to stream.
+ *
+ * If vector is larger than 1000 nothing is printed.
  *
  * @param stream  Prints only, if stream is active
  * @param name    Specific string to print (usually name of the vector)
@@ -452,14 +509,19 @@ void printVector_gb(enum LOG_STREAM stream, char name[], double* a, int n, doubl
   // This only works for number of states less than 10!
   // For large arrays, this is not a good output format!
   char row_to_print[40960];
-  sprintf(row_to_print, "%s(%8g) =\t", name, time);
-  for (int i=0;i<n;i++)
-    sprintf(row_to_print, "%s %18.12g", row_to_print, a[i]);
+  unsigned int bufSize = 40960;
+  unsigned int ct;
+  ct += snprintf(row_to_print, bufSize, "%s(%8g) =\t", name, time);
+  for (int i=0;i<n;i++) {
+    ct += snprintf(row_to_print+ct, bufSize-ct, "%18.12g", a[i]);
+  }
   infoStreamPrint(stream, 0, "%s", row_to_print);
 }
 
 /**
- * @brief Prints an integer vector
+ * @brief Prints an integer vector to stream.
+ *
+ * If vector is larger than 1000 nothing is printed.
  *
  * @param name    Specific string to print (usually name of the vector)
  * @param a       Integer vector to print
@@ -472,33 +534,18 @@ void printIntVector_gb(enum LOG_STREAM stream, char name[], int* a, int n, doubl
   if (!ACTIVE_STREAM(stream) || n>1000) return;
 
   char row_to_print[40960];
-  sprintf(row_to_print, "%s(%8g) =\t", name, time);
+  unsigned int bufSize = 40960;
+  unsigned int ct;
+  ct = snprintf(row_to_print, bufSize, "%s(%8g) =\t", name, time);
   for (int i=0;i<n;i++)
-    sprintf(row_to_print, "%s %d", row_to_print, a[i]);
+    ct += snprintf(row_to_print+ct, bufSize-ct, "%d", a[i]);
   infoStreamPrint(stream, 0, "%s", row_to_print);
 }
 
 /**
- * @brief Prints a square matrix
+ * @brief Prints selected vector components given by an index vector.
  *
- * @param name    Specific string to print (usually name of the matrix)
- * @param a       Matrix to print
- * @param n       number of columns and rows
- * @param time    Time value
- */
-void printMatrix_gb(char name[], double* a, int n, double time) {
-  printf("\n%s at time: %g: \n ", name, time);
-  for (int i=0;i<n;i++)
-  {
-    for (int j=0;j<n;j++)
-      printf("%6g ", a[i*n + j]);
-    printf("\n");
-  }
-  printf("\n");
-}
-
-/**
- * @brief Prints selected vector components given by an index vector
+ * If more than 1000 elements should be printed do nothing.
  *
  * @param name    Specific string to print (usually name of the vector)
  * @param a       Vector to print
@@ -515,9 +562,11 @@ void printVector_gbf(enum LOG_STREAM stream, char name[], double* a, int n, doub
   // This only works for number of states less than 10!
   // For large arrays, this is not a good output format!
   char row_to_print[40960];
-  sprintf(row_to_print, "%s(%8g) =\t", name, time);
+  unsigned int bufSize = 40960;
+  unsigned int ct;
+  ct = snprintf(row_to_print, bufSize, "%s(%8g) =\t", name, time);
   for (int i=0;i<nIndx;i++)
-    sprintf(row_to_print, "%s %16.12g", row_to_print, a[indx[i]]);
+    ct += snprintf(row_to_print+ct, bufSize-ct, "%16.12g", a[indx[i]]);
   infoStreamPrint(stream, 0, "%s", row_to_print);
 }
 
@@ -565,39 +614,55 @@ void printSparseJacobianLocal(ANALYTIC_JACOBIAN* jacobian, const char* name) {
 /**
  * @brief Write information on the active fast states on file (activity diagram)
  *
- * @param gbData  Pointer to generik GBODE data struct.
- * @param event   If an event has happend, write zeros else ones
- * @param time    Actual time of reporting
+ * @param gbData       Pointer to generic GBODE data struct.
+ * @param event        If an event has happened, write zeros else ones
+ * @param time         Actual time of reporting
+ * @param rejectedType Type of rejection
+ *                     0  <= no rejection
+ *                     1  <= error of slow states greater than the tolerance
+ *                     2  <= interpolation error is too large
+ *                     3  <= rejected because solving the NLS failed
+ *                    -1  <= step is preliminary accepted but needs refinement
  */
-void dumpFastStates_gb(DATA_GBODE* gbData, modelica_boolean event, double time) {
-    char fastStates_row[4096];
-    sprintf(fastStates_row, "%15.10g %15.10g %15.10g %15.10g", time, gbData->err_slow, gbData->err_int, gbData->err_fast);
-    for (int i = 0; i < gbData->nStates; i++) {
-      if (event)
-        sprintf(fastStates_row, "%s 0", fastStates_row);
-      else
-        sprintf(fastStates_row, "%s 1", fastStates_row);
-    }
-    fprintf(gbData->gbfData->fastStatesDebugFile, "%s\n", fastStates_row);
+void dumpFastStates_gb(DATA_GBODE* gbData, modelica_boolean event, double time, int rejectedType) {
+  char fastStates_row[4096];
+  unsigned int bufSize = 4096;
+  unsigned int ct;
+  ct = snprintf(fastStates_row, bufSize, "%15.10g %2d %15.10g %15.10g %15.10g", time, rejectedType, gbData->err_slow, gbData->err_int, gbData->err_fast);
+  for (int i = 0; i < gbData->nStates; i++) {
+    if (event)
+      ct += snprintf(fastStates_row+ct, bufSize-ct, " 0");
+    else
+      ct += snprintf(fastStates_row+ct, bufSize-ct, " 1");
+  }
+  fprintf(gbData->gbfData->fastStatesDebugFile, "%s\n", fastStates_row);
 }
 
 /**
  * @brief Write information on the active fast states on file (activity diagram)
  *
- * @param gbData  Pointer to generik GBODE data struct.
+ * @param gbData  Pointer to generic GBODE data struct.
  * @param time    Actual time of reporting
+ * @param rejectedType Type of rejection
+ *                     0  <= no rejection
+ *                     1  <= error of slow states greater than the tolerance
+ *                     2  <= interpolation error is too large
+ *                     3  <= rejected because solving the NLS failed
+ *                    -1  <= step is preliminary accepted but needs refinement
  */
-void dumpFastStates_gbf(DATA_GBODE* gbData, double time) {
+void dumpFastStates_gbf(DATA_GBODE* gbData, double time, int rejectedType) {
   char fastStates_row[4096];
+  unsigned int bufSize = 40960;
+  unsigned int ct;
   int i, ii;
-  sprintf(fastStates_row, "%15.10g %15.10g %15.10g %15.10g", time, gbData->err_slow, gbData->err_int, gbData->err_fast);
+  ct = snprintf(fastStates_row, bufSize, "%15.10g %2d %15.10g %15.10g %15.10g", time, rejectedType, gbData->err_slow, gbData->err_int, gbData->err_fast);
   for (i = 0, ii = 0; i < gbData->nStates;) {
     if (i == gbData->fastStatesIdx[ii]) {
-      sprintf(fastStates_row, "%s 1", fastStates_row);
+      ct += snprintf(fastStates_row+ct, bufSize-ct, " 1");
       i++;
-      ii++;
+      if (ii < gbData->nFastStates-1) ii++;
     } else {
-      sprintf(fastStates_row, "%s 0", fastStates_row);
+      ct += snprintf(fastStates_row+ct, bufSize-ct, " 0");
       i++;
     }
   }
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_util.h b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_util.h
index 314a19df90f..1b3f69bff72 100644
--- a/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_util.h
+++ b/OMCompiler/SimulationRuntime/c/simulation/solver/gbode_util.h
@@ -64,13 +64,12 @@ void projVector_gbf(double* a, double* b, int nIndx, int* indx);
 // Debug functions for the development of gbode
 void printVector_gb(enum LOG_STREAM stream, char name[], double* a, int n, double time);
 void printIntVector_gb(enum LOG_STREAM stream, char name[], int* a, int n, double time);
-void printMatrix_gb(char name[], double* a, int n, double time);
 void printVector_gbf(enum LOG_STREAM stream, char name[], double* a, int n, double time, int nIndx, int* indx);
 void printSparseJacobianLocal(ANALYTIC_JACOBIAN* jacobian, const char* name);
 
 void debugRingBuffer(enum LOG_STREAM stream, double* x, double* k, int nStates, BUTCHER_TABLEAU* tableau, double time, double stepSize);
-void dumpFastStates_gb(DATA_GBODE *gbData, modelica_boolean event, double time);
-void dumpFastStates_gbf(DATA_GBODE *gbData, double time);
+void dumpFastStates_gb(DATA_GBODE *gbData, modelica_boolean event, double time, int rejectedType);
+void dumpFastStates_gbf(DATA_GBODE *gbData, double time, int rejectedType);
 
 modelica_boolean checkFastStatesChange(DATA_GBODE* gbData);
 
diff --git a/OMCompiler/SimulationRuntime/c/simulation/solver/solver_main.c b/OMCompiler/SimulationRuntime/c/simulation/solver/solver_main.c
index f2104c75b82..e09c3170e33 100644
--- a/OMCompiler/SimulationRuntime/c/simulation/solver/solver_main.c
+++ b/OMCompiler/SimulationRuntime/c/simulation/solver/solver_main.c
@@ -192,7 +192,7 @@ int solver_main_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverIn
       data->simulationInfo->solverSteps = solverInfo->solverStats[0] + solverInfo->solverStatsTmp[0];
     return retVal;
   default:
-    throwStreamPrint(threadData, "Unhandles case in solver_main_step.");
+    throwStreamPrint(threadData, "Unhandled case in solver_main_step.");
   }
   TRACE_POP
   return 1;
diff --git a/testsuite/simulation/modelica/solver/gbode/IRKGaussian_01.mos b/testsuite/simulation/modelica/solver/gbode/IRKGaussian_01.mos
index 043aeba0db5..42a1a7ed045 100644
--- a/testsuite/simulation/modelica/solver/gbode/IRKGaussian_01.mos
+++ b/testsuite/simulation/modelica/solver/gbode/IRKGaussian_01.mos
@@ -36,7 +36,7 @@ rkMethods := {"gauss2",
 
 nlsMethods := {"kinsol"};
 
-errCtrls := {"0"};
+errCtrls := {"const"};
 
 setCommandLineOptions("--generateSymbolicJacobian"); getErrorString();
 
@@ -80,7 +80,7 @@ end for;
 // ""
 // {"gauss2","gauss3","gauss4","gauss5","gauss6","radauIA2","radauIA3","radauIA4","radauIIA2","radauIIA3","radauIIA4","lobattoIIIA3","lobattoIIIA4","lobattoIIIB3","lobattoIIIB4","lobattoIIIC3","lobattoIIIC4"}
 // {"kinsol"}
-// {"0"}
+// {"const"}
 // true
 // ""
 // {"SlowFastDynamics","SlowFastDynamics_init.xml"}
diff --git a/testsuite/simulation/modelica/solver/gbode/IRK_01.mos b/testsuite/simulation/modelica/solver/gbode/IRK_01.mos
index d8a41341052..ebc753a523f 100644
--- a/testsuite/simulation/modelica/solver/gbode/IRK_01.mos
+++ b/testsuite/simulation/modelica/solver/gbode/IRK_01.mos
@@ -24,7 +24,7 @@ rkMethods := {"impl_euler",
 
 nlsMethods := {"newton", "kinsol"};
 
-errCtrls := {"0", "1"};
+errCtrls := {"const", "i"};
 
 setCommandLineOptions("--generateSymbolicJacobian"); getErrorString();
 
@@ -68,7 +68,7 @@ end for;
 // ""
 // {"impl_euler","sdirk2","sdirk3","esdirk2","esdirk3"}
 // {"newton","kinsol"}
-// {"0","1"}
+// {"const","i"}
 // true
 // ""
 // {"SlowFastDynamics","SlowFastDynamics_init.xml"}
diff --git a/testsuite/simulation/modelica/solver/gbode/RK_01.mos b/testsuite/simulation/modelica/solver/gbode/RK_01.mos
index 653a12f6ba5..327f53ce245 100644
--- a/testsuite/simulation/modelica/solver/gbode/RK_01.mos
+++ b/testsuite/simulation/modelica/solver/gbode/RK_01.mos
@@ -26,7 +26,7 @@ rkMethods := {"expl_euler",
              "rk1012",
              "rk1214"};
 
-errCtrls := {"0", "1"};
+errCtrls := {"const", "i"};
 
 setCommandLineOptions("--generateSymbolicJacobian"); getErrorString();
 
@@ -67,7 +67,7 @@ end for;
 // true
 // ""
 // {"expl_euler","dopri45","merson","fehlberg12","fehlberg45","fehlberg78","rk810","rk1012","rk1214"}
-// {"0","1"}
+// {"const","i"}
 // true
 // ""
 // {"SlowFastDynamics","SlowFastDynamics_init.xml"}
diff --git a/testsuite/simulation/modelica/solver/gbode/multiRate_01.mos b/testsuite/simulation/modelica/solver/gbode/multiRate_01.mos
index f78293a8b73..76270273fab 100644
--- a/testsuite/simulation/modelica/solver/gbode/multiRate_01.mos
+++ b/testsuite/simulation/modelica/solver/gbode/multiRate_01.mos
@@ -21,10 +21,11 @@ rkMethods := {"expl_euler",
               "sdirk3",
               "esdirk2",
               "esdirk3",
+              "esdirk4",
               "merson",
-              "dopri45"};
+              "dopri45 -gbint=dense_output_errctrl"};
 
-nlsMethods := {"newton", "kinsol"};
+nlsMethods := {"newton"};
 
 setCommandLineOptions("--generateSymbolicJacobian"); getErrorString();
 
@@ -63,8 +64,8 @@ end for;
 // Result:
 // true
 // ""
-// {"expl_euler","impl_euler","sdirk3","esdirk2","esdirk3","merson","dopri45"}
-// {"newton","kinsol"}
+// {"expl_euler","impl_euler","sdirk3","esdirk2","esdirk3","esdirk4","merson","dopri45 -gbint=dense_output_errctrl"}
+// {"newton"}
 // true
 // ""
 // {"SlowFastDynamics","SlowFastDynamics_init.xml"}
@@ -88,16 +89,6 @@ end for;
 // Failed vars:
 //
 // --------------------------------------------------------
-// Running RK expl_euler with NLS kinsol:
-// stdout            | warning | Start or stop time was overwritten, but no new integrator step size was provided.
-// |                 | info    | | Re-calculating step size for 500 intervals.
-// |                 | |       | | Add `stepSize=<value>` to `-override=` or override file to silence this warning.
-// LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
-// LOG_SUCCESS       | info    | The simulation finished successfully.
-//
-// Failed vars:
-//
-// --------------------------------------------------------
 // Running RK impl_euler with NLS newton:
 // stdout            | warning | Start or stop time was overwritten, but no new integrator step size was provided.
 // |                 | info    | | Re-calculating step size for 500 intervals.
@@ -108,16 +99,6 @@ end for;
 // Failed vars:
 //
 // --------------------------------------------------------
-// Running RK impl_euler with NLS kinsol:
-// stdout            | warning | Start or stop time was overwritten, but no new integrator step size was provided.
-// |                 | info    | | Re-calculating step size for 500 intervals.
-// |                 | |       | | Add `stepSize=<value>` to `-override=` or override file to silence this warning.
-// LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
-// LOG_SUCCESS       | info    | The simulation finished successfully.
-//
-// Failed vars:
-//
-// --------------------------------------------------------
 // Running RK sdirk3 with NLS newton:
 // stdout            | warning | Start or stop time was overwritten, but no new integrator step size was provided.
 // |                 | info    | | Re-calculating step size for 500 intervals.
@@ -128,16 +109,6 @@ end for;
 // Failed vars:
 //
 // --------------------------------------------------------
-// Running RK sdirk3 with NLS kinsol:
-// stdout            | warning | Start or stop time was overwritten, but no new integrator step size was provided.
-// |                 | info    | | Re-calculating step size for 500 intervals.
-// |                 | |       | | Add `stepSize=<value>` to `-override=` or override file to silence this warning.
-// LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
-// LOG_SUCCESS       | info    | The simulation finished successfully.
-//
-// Failed vars:
-//
-// --------------------------------------------------------
 // Running RK esdirk2 with NLS newton:
 // stdout            | warning | Start or stop time was overwritten, but no new integrator step size was provided.
 // |                 | info    | | Re-calculating step size for 500 intervals.
@@ -148,16 +119,6 @@ end for;
 // Failed vars:
 //
 // --------------------------------------------------------
-// Running RK esdirk2 with NLS kinsol:
-// stdout            | warning | Start or stop time was overwritten, but no new integrator step size was provided.
-// |                 | info    | | Re-calculating step size for 500 intervals.
-// |                 | |       | | Add `stepSize=<value>` to `-override=` or override file to silence this warning.
-// LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
-// LOG_SUCCESS       | info    | The simulation finished successfully.
-//
-// Failed vars:
-//
-// --------------------------------------------------------
 // Running RK esdirk3 with NLS newton:
 // stdout            | warning | Start or stop time was overwritten, but no new integrator step size was provided.
 // |                 | info    | | Re-calculating step size for 500 intervals.
@@ -168,7 +129,7 @@ end for;
 // Failed vars:
 //
 // --------------------------------------------------------
-// Running RK esdirk3 with NLS kinsol:
+// Running RK esdirk4 with NLS newton:
 // stdout            | warning | Start or stop time was overwritten, but no new integrator step size was provided.
 // |                 | info    | | Re-calculating step size for 500 intervals.
 // |                 | |       | | Add `stepSize=<value>` to `-override=` or override file to silence this warning.
@@ -188,27 +149,7 @@ end for;
 // Failed vars:
 //
 // --------------------------------------------------------
-// Running RK merson with NLS kinsol:
-// stdout            | warning | Start or stop time was overwritten, but no new integrator step size was provided.
-// |                 | info    | | Re-calculating step size for 500 intervals.
-// |                 | |       | | Add `stepSize=<value>` to `-override=` or override file to silence this warning.
-// LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
-// LOG_SUCCESS       | info    | The simulation finished successfully.
-//
-// Failed vars:
-//
-// --------------------------------------------------------
-// Running RK dopri45 with NLS newton:
-// stdout            | warning | Start or stop time was overwritten, but no new integrator step size was provided.
-// |                 | info    | | Re-calculating step size for 500 intervals.
-// |                 | |       | | Add `stepSize=<value>` to `-override=` or override file to silence this warning.
-// LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
-// LOG_SUCCESS       | info    | The simulation finished successfully.
-//
-// Failed vars:
-//
-// --------------------------------------------------------
-// Running RK dopri45 with NLS kinsol:
+// Running RK dopri45 -gbint=dense_output_errctrl with NLS newton:
 // stdout            | warning | Start or stop time was overwritten, but no new integrator step size was provided.
 // |                 | info    | | Re-calculating step size for 500 intervals.
 // |                 | |       | | Add `stepSize=<value>` to `-override=` or override file to silence this warning.