diff --git a/src/OMSimulatorLib/AlgLoop.cpp b/src/OMSimulatorLib/AlgLoop.cpp
index 033f98696..26a800d96 100644
--- a/src/OMSimulatorLib/AlgLoop.cpp
+++ b/src/OMSimulatorLib/AlgLoop.cpp
@@ -29,6 +29,10 @@
  *
  */
 
+#include <sstream>
+#include <cmath>
+#include <algorithm>
+
 #include "AlgLoop.h"
 
 #include "Component.h"
@@ -36,7 +40,13 @@
 #include "Flags.h"
 #include "System.h"
 
-#include <sstream>
+#ifdef min
+#undef min
+#endif
+
+#ifdef max
+#undef max
+#endif
 
 /**
  * @brief Check flag returned by KINSOL function and log error
@@ -292,7 +302,7 @@ oms::KinsolSolver::~KinsolSolver()
  * @param absoluteTolerance     Tolerance used for solving the loop
  * @return oms::KinsolSolver*   Retruns pointer to KinsolSolver object
  */
-oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance, const bool useDirectionalDerivative)
+oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance, double relativeTolerance, const bool useDirectionalDerivative)
 {
   int flag;
   int printLevel;
@@ -301,6 +311,7 @@ oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, cons
   logDebug("Create new KinsolSolver object for algebraic loop number " + std::to_string(algLoopNum));
 
   kinsolSolver->size = size;
+  kinsolSolver->firstSolution = true;
 
   /* Allocate memory */
   kinsolSolver->initialGuess = N_VNew_Serial(kinsolSolver->size);
@@ -364,10 +375,12 @@ oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, cons
   if (!checkFlag(flag, "KINSetJacFn")) return NULL;
 
   /* Set function-norm stopping tolerance */
-  kinsolSolver->fnormtol = absoluteTolerance;
+  kinsolSolver->fnormtol = 0.01 * absoluteTolerance;
   flag = KINSetFuncNormTol(kinsolSolver->kinsolMemory, kinsolSolver->fnormtol);
   if (!checkFlag(flag, "KINSetFuncNormTol")) return NULL;
 
+  kinsolSolver->freltol = relativeTolerance * 0.01;
+
   /* Set scaled-step stopping tolerance */
   flag = KINSetScaledStepTol(kinsolSolver->kinsolMemory, 0.0);
   if (!checkFlag(flag, "KINSetScaledStepTol")) return NULL;
@@ -400,7 +413,7 @@ oms::KinsolSolver* oms::KinsolSolver::NewKinsolSolver(const int algLoopNum, cons
  *                            `oms_status_warning` if solving was computed, but solution is not within tolerance
  *                            and `oms_status_error` if an error occured.
  */
-oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& graph)
+oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& graph, double tolerance /*= 0.0*/)
 {
   /* Update user data */
   KINSOL_USER_DATA* kinsolUserData = (KINSOL_USER_DATA*) user_data;
@@ -413,7 +426,7 @@ oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& gra
   int flag;
   double fNormValue;
 
-  logDebug("Solving system " + std::to_string(kinsolUserData->algLoopNumber));
+  double tol = tolerance != 0.0 ? tolerance : fnormtol;
 
   if (SCC.connections.size() != size)
   {
@@ -432,9 +445,29 @@ oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& gra
     }
   }
 
+  /* Apply relative tolerance to magnitude of initial guess */
+  fNormValue = std::sqrt(N_VDotProd_Serial(initialGuess, initialGuess));
+  tol = std::max(tol, fNormValue * freltol * tolerance / fnormtol);
+
+  /* Check residual for initial guess */
+  nlsKinsolResiduals(initialGuess, fTmp, user_data);
+  fNormValue = std::sqrt(N_VDotProd_Serial(fTmp, fTmp));
+  if (fNormValue > 0.0)
+    tol = std::min(fNormValue, tol); // We already know this is achievable, but let's get in at least one iteration
+  else
+    return oms_status_ok;
+
+  if (!firstSolution) {
+    // Predict a better initial guess
+    SUNLinSolSolve(linSol, J, y, fTmp, tol);
+    N_VLinearSum(1.0, initialGuess, -1.0, y, initialGuess);
+  }
+
   /* u and f scaling */
   // TODO: Add scaling that is not only constant ones
 
+  KINSetFuncNormTol(kinsolMemory, tol);
+
   /* Solve algebraic loop with KINSol() */
   flag = KINSol(kinsolMemory,   /* KINSol memory block */
                 initialGuess,   /* initial guess on input; solution vector */
@@ -444,15 +477,19 @@ oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& gra
   if (!checkFlag(flag, "KINSol")) return oms_status_error;
 
   /* Check solution */
-  flag = nlsKinsolResiduals(initialGuess, fTmp, user_data);
-  fNormValue = N_VWL2Norm(fTmp, fTmp);
-  if ( fNormValue > fnormtol )
+  KINGetFuncNorm(kinsolMemory, &fNormValue);
+  if ( fNormValue > tol )
   {
     logWarning("Solution of algebraic loop " + std::to_string(((KINSOL_USER_DATA *)user_data)->algLoopNumber) + "not within precission given by fnormtol: " + std::to_string(fnormtol));
     logDebug("2-norm of residual of solution: " + std::to_string(fNormValue));
     return oms_status_warning;
   }
 
+  if (flag == KIN_SUCCESS) {
+    firstSolution = false;
+    KINSetNoInitSetup(kinsolMemory, SUNTRUE);
+  }
+
   logDebug("Solved system " + std::to_string(kinsolUserData->algLoopNumber) + " successfully");
   return oms_status_ok;
 }
@@ -467,7 +504,7 @@ oms_status_enu_t oms::KinsolSolver::kinsolSolve(System& syst, DirectedGraph& gra
  * @param scc     Strong Connected Componten, a vector of connected
  * @param number
  */
-oms::AlgLoop::AlgLoop(oms_alg_solver_enu_t method, double absTol, scc_t scc, const int number, const bool useDirectionalDerivative): absoluteTolerance(absTol), SCC(scc), systNumber(number)
+oms::AlgLoop::AlgLoop(oms_alg_solver_enu_t method, double absTol, double relTol, scc_t scc, const int number, const bool useDirectionalDerivative): absoluteTolerance(absTol), SCC(scc), systNumber(number)
 {
   switch (method)
   {
@@ -482,7 +519,7 @@ oms::AlgLoop::AlgLoop(oms_alg_solver_enu_t method, double absTol, scc_t scc, con
 
   if (method == oms_alg_solver_kinsol)
   {
-    kinsolData = KinsolSolver::NewKinsolSolver(systNumber, SCC.connections.size(), absoluteTolerance, useDirectionalDerivative);
+    kinsolData = KinsolSolver::NewKinsolSolver(systNumber, SCC.connections.size(), absoluteTolerance, relTol, useDirectionalDerivative);
     if (kinsolData==NULL)
     {
       logError("NewKinsolSolver() failed. Aborting!");
@@ -501,7 +538,7 @@ oms::AlgLoop::AlgLoop(oms_alg_solver_enu_t method, double absTol, scc_t scc, con
  * @param graph               Reference to directed graph
  * @return oms_status_enu_t   Returns `oms_status_ok` on success
  */
-oms_status_enu_t oms::AlgLoop::solveAlgLoop(System& syst, DirectedGraph& graph)
+oms_status_enu_t oms::AlgLoop::solveAlgLoop(System& syst, DirectedGraph& graph, double tolerance)
 {
   logDebug("Solving algebraic loop formed by connections\n" + dumpLoopVars(graph));
   logDebug("Using solver " + getAlgSolverName());
@@ -509,9 +546,9 @@ oms_status_enu_t oms::AlgLoop::solveAlgLoop(System& syst, DirectedGraph& graph)
   switch (algSolverMethod)
   {
   case oms_alg_solver_fixedpoint:
-    return fixPointIteration(syst, graph);
+    return fixPointIteration(syst, graph, tolerance);
   case oms_alg_solver_kinsol:
-    return kinsolData->kinsolSolve(syst, graph);
+    return kinsolData->kinsolSolve(syst, graph, tolerance);
   default:
     logError("Invalid algebraic solver method!");
     return oms_status_error;
@@ -525,7 +562,7 @@ oms_status_enu_t oms::AlgLoop::solveAlgLoop(System& syst, DirectedGraph& graph)
  * @param graph
  * @return oms_status_enu_t
  */
-oms_status_enu_t oms::AlgLoop::fixPointIteration(System& syst, DirectedGraph& graph)
+oms_status_enu_t oms::AlgLoop::fixPointIteration(System& syst, DirectedGraph& graph, double tolerance)
 {
   const int size = SCC.connections.size();
   const int maxIterations = Flags::MaxLoopIteration();
@@ -533,6 +570,9 @@ oms_status_enu_t oms::AlgLoop::fixPointIteration(System& syst, DirectedGraph& gr
   double maxRes;
   double *res = new double[size]();
 
+  if (tolerance == 0.0 || tolerance > absoluteTolerance)
+    tolerance = absoluteTolerance;
+
   do
   {
     std::stringstream ss;
@@ -641,7 +681,7 @@ oms_status_enu_t oms::AlgLoop::fixPointIteration(System& syst, DirectedGraph& gr
       logInfo(ss.str());
     }
 
-  } while(maxRes > absoluteTolerance && it < maxIterations);
+  } while(maxRes > tolerance && it < maxIterations);
 
   delete[] res;
 
diff --git a/src/OMSimulatorLib/AlgLoop.h b/src/OMSimulatorLib/AlgLoop.h
index f25d945e8..fa9b47d73 100644
--- a/src/OMSimulatorLib/AlgLoop.h
+++ b/src/OMSimulatorLib/AlgLoop.h
@@ -57,12 +57,13 @@ namespace oms
   {
   public:
     ~KinsolSolver();
-    static KinsolSolver* NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance, const bool useDirectionalDerivative);
-    oms_status_enu_t kinsolSolve(System& syst, DirectedGraph& graph);
+    static KinsolSolver* NewKinsolSolver(const int algLoopNum, const unsigned int size, double absoluteTolerance, double relativeTolerance, const bool useDirectionalDerivative);
+    oms_status_enu_t kinsolSolve(System& syst, DirectedGraph& graph, double tolerance = 0.0);
 
   private:
     /* tolerances */
     double fnormtol;        /* function tolerance */
+    double freltol;         /* relative function tolerance */
 
     /* work arrays */
     N_Vector initialGuess;
@@ -80,6 +81,8 @@ namespace oms
     N_Vector y;             /* Template for cloning vectors needed inside linear solver */
     SUNMatrix J;            /* (Non-)Sparse matrix template for cloning matrices needed within linear solver */
 
+    bool firstSolution;
+
     /* member function */
     static int nlsKinsolJac(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2);
     static int nlsKinsolResiduals(N_Vector u, N_Vector fval, void *user_data);
@@ -90,16 +93,16 @@ namespace oms
   class AlgLoop
   {
   public:
-    AlgLoop(oms_alg_solver_enu_t method, double absTol, scc_t SCC, const int systNumber, const bool useDirectionalDerivative);
+    AlgLoop(oms_alg_solver_enu_t method, double absTol, double relTol, scc_t SCC, const int systNumber, const bool useDirectionalDerivative);
 
     scc_t getSCC() {return SCC;}
-    oms_status_enu_t solveAlgLoop(System& syst, DirectedGraph& graph);
+    oms_status_enu_t solveAlgLoop(System& syst, DirectedGraph& graph, double tolerance);
     std::string getAlgSolverName();
     std::string dumpLoopVars(DirectedGraph& graph);
 
   private:
     oms_alg_solver_enu_t algSolverMethod;
-    oms_status_enu_t fixPointIteration(System& syst, DirectedGraph& graph);
+    oms_status_enu_t fixPointIteration(System& syst, DirectedGraph& graph, double tolerance);
 
     KinsolSolver* kinsolData;
 
diff --git a/src/OMSimulatorLib/System.cpp b/src/OMSimulatorLib/System.cpp
index 0c99305a3..8b5a24057 100644
--- a/src/OMSimulatorLib/System.cpp
+++ b/src/OMSimulatorLib/System.cpp
@@ -3066,7 +3066,7 @@ oms_status_enu_t oms::System::addAlgLoop(scc_t SCC, const int algLoopNum, Direct
     loopsNeedUpdate = false;
   }
 
-  algLoops.push_back( AlgLoop(Flags::AlgLoopSolver(), absoluteTolerance, SCC, algLoopNum, supportsDirectionalDerivatives));
+  algLoops.push_back( AlgLoop(Flags::AlgLoopSolver(), absoluteTolerance, relativeTolerance, SCC, algLoopNum, supportsDirectionalDerivatives));
 
   return oms_status_ok;
 }
@@ -3110,9 +3110,9 @@ oms_status_enu_t oms::System::updateAlgebraicLoops(const std::vector<scc_t>& sor
 }
 
 
-oms_status_enu_t oms::System::solveAlgLoop(DirectedGraph& graph, int loopNumber)
+oms_status_enu_t oms::System::solveAlgLoop(DirectedGraph& graph, int loopNumber, double tolerance)
 {
-  return algLoops[loopNumber].solveAlgLoop(*this, graph);
+  return algLoops[loopNumber].solveAlgLoop(*this, graph, tolerance);
 }
 
 oms_status_enu_t oms::System::rename(const oms::ComRef& newCref)
diff --git a/src/OMSimulatorLib/System.h b/src/OMSimulatorLib/System.h
index bb1d3562e..0775f30b2 100644
--- a/src/OMSimulatorLib/System.h
+++ b/src/OMSimulatorLib/System.h
@@ -188,7 +188,7 @@ namespace oms
     AlgLoop* getAlgLoop(const int systemNumber);
     oms_status_enu_t addAlgLoop(scc_t SCC, const int algLoopNum, DirectedGraph& graph, bool supportsDirectionalDerivatives);
     oms_status_enu_t updateAlgebraicLoops(const std::vector< scc_t >& sortedConnections, DirectedGraph& graph);
-    oms_status_enu_t solveAlgLoop(DirectedGraph& graph, int loopNumber);
+    oms_status_enu_t solveAlgLoop(DirectedGraph& graph, int loopNumber, double tolerance = 0.0);
 
     bool useThreadPool();
     ctpl::thread_pool& getThreadPool();
diff --git a/src/OMSimulatorLib/SystemSC.cpp b/src/OMSimulatorLib/SystemSC.cpp
index 157ff0f40..da1ea0f27 100644
--- a/src/OMSimulatorLib/SystemSC.cpp
+++ b/src/OMSimulatorLib/SystemSC.cpp
@@ -44,27 +44,55 @@
 
 int oms::cvode_rhs(realtype t, N_Vector y, N_Vector ydot, void* user_data)
 {
+  static const double U = sqrt(DBL_EPSILON);
+
   SystemSC* system = (SystemSC*)user_data;
   oms_status_enu_t status;
   fmi2Status fmistatus;
 
+  // smallest step made by the difference quotient Jacobian evaluation
+  double sigmaMin = 0.0;
+  double fnorm = 0.0;
+
+  double h = 0.0;
+  CVodeGetCurrentStep(system->solverData.cvode.mem, &h);
+
+  size_t maxStates = system->nStates.empty() ? 0 : *std::max_element(system->nStates.begin(), system->nStates.end());
+  std::unique_ptr<double> vs(new double[maxStates]);
+  double* values = vs.get();
+
   // update states in FMUs
-  for (size_t i=0, j=0; i < system->fmus.size(); ++i)
+  size_t j = 0;
+  for (size_t i=0; i < system->fmus.size(); ++i)
   {
     system->fmus[i]->setTime(t);
 
     if (0 == system->nStates[i])
       continue;
 
-    for (size_t k = 0; k < system->nStates[i]; k++, j++)
-      system->states[i][k] = NV_Ith_S(y, j);
+    for (size_t k = 0; k < system->nStates[i]; k++, j++) {
+      values[k] = NV_Ith_S(y, j);
+
+      double w = NV_Ith_S(system->solverData.cvode.ewt, j);
+      double f = system->states_der[i][k] * w;
+      fnorm += f * f;
+
+      // smallest possible difference quotient step size used by CVODE
+      double s = 1.0 / w;
+      if (sigmaMin == 0.0 || s < sigmaMin)
+        sigmaMin = s;
+    }
 
     // set states
-    status = system->fmus[i]->setContinuousStates(system->states[i]);
+    status = system->fmus[i]->setContinuousStates(values);
     if (oms_status_ok != status) return status;
   }
 
-  system->updateInputs(system->eventGraph);
+  if (fnorm > 0.0 && h > 0.0)
+    sigmaMin *= 1000.0 * DBL_EPSILON * h * j * sqrt(fnorm);
+
+  // Solve loops to a sufficient precision for Jacobian calculation
+  system->updateInputs(system->eventGraph, 8 * sigmaMin);
 
   // get state derivatives
   for (size_t j=0, k=0; j < system->fmus.size(); ++j)
@@ -73,11 +101,11 @@ int oms::cvode_rhs(realtype t, N_Vector y, N_Vector ydot, void* user_data)
       continue;
 
     // get state derivatives
-    status = system->fmus[j]->getDerivatives(system->states_der[j]);
+    status = system->fmus[j]->getDerivatives(values);
     if (oms_status_ok != status) return status;
 
     for (size_t i=0; i < system->nStates[j]; ++i, ++k)
-      NV_Ith_S(ydot, k) = system->states_der[j][i];
+      NV_Ith_S(ydot, k) = values[i];
   }
 
   return 0;
@@ -347,9 +375,15 @@ oms_status_enu_t oms::SystemSC::initialize()
 
     solverData.cvode.abstol = N_VNew_Serial(static_cast<long>(n_states));
     if (!solverData.cvode.abstol) logError("SUNDIALS_ERROR: N_VNew_Serial() failed - returned NULL pointer");
-    for (size_t j=0, k=0; j < fmus.size(); ++j)
-      for (size_t i=0; i < nStates[j]; ++i, ++k)
-        NV_Ith_S(solverData.cvode.abstol, k) = 0.01*absoluteTolerance*states_nominal[j][i];
+    solverData.cvode.ewt = N_VNew_Serial(static_cast<long>(n_states));
+    if (!solverData.cvode.ewt) logError("SUNDIALS_ERROR: N_VClone() failed - returned NULL pointer");
+    for (size_t j=0, k=0; j < fmus.size(); ++j) {
+      for (size_t i=0; i < nStates[j]; ++i, ++k) {
+        double abstol = absoluteTolerance * states_nominal[j][i];
+        NV_Ith_S(solverData.cvode.abstol, k) = abstol;
+        NV_Ith_S(solverData.cvode.ewt, k) = 1.0 / (relativeTolerance * fabs(states[j][i]) + abstol);
+      }
+    }
     //N_VPrint_Serial(solverData.cvode.abstol);
 
     // Call CVodeCreate to create the solver memory and specify the
@@ -455,6 +489,7 @@ oms_status_enu_t oms::SystemSC::terminate()
     SUNLinSolFree(solverData.cvode.linSol);
     N_VDestroy_Serial(solverData.cvode.y);
     N_VDestroy_Serial(solverData.cvode.abstol);
+    N_VDestroy_Serial(solverData.cvode.ewt);
     CVodeFree(&(solverData.cvode.mem));
     solverData.cvode.mem = NULL;
   }
@@ -534,23 +569,23 @@ oms_status_enu_t oms::SystemSC::doStep()
   switch(solverMethod)
   {
     case oms_solver_sc_explicit_euler:
-      return doStepEuler();
+      return doStepEuler(time + maximumStepSize);
 
     case oms_solver_sc_cvode:
-      return doStepCVODE();
+      return doStepCVODE(time + maximumStepSize);
 
     default:
       return logError_InternalError;
   }
 }
 
-oms_status_enu_t oms::SystemSC::doStepEuler()
+oms_status_enu_t oms::SystemSC::doStepEuler(double stopTime)
 {
   fmi2Status fmistatus;
   oms_status_enu_t status;
 
   // Step 1: Initialize state variables and time
-  const fmi2Real end_time = std::min(time + maximumStepSize, getModel().getStopTime());
+  const fmi2Real end_time = std::min(time + maximumStepSize, stopTime);
   const fmi2Real event_time_tolerance = 1e-4;
 
   logDebug("doStepEuler: " + std::to_string(time) + " -> " + std::to_string(end_time));
@@ -613,6 +648,10 @@ oms_status_enu_t oms::SystemSC::doStepEuler()
     step_size_adjustment *= 0.5; // reduce the step size in each iteration
 
     // a. Evaluate derivatives for each FMU
+    // set time
+    for (const auto& component : getComponents())
+      component.second->setTime(time);
+      
     const fmi2Real step_size = event_time - time;  // Substep size, do one step from current time to the event
     logDebug("step_size: " + std::to_string(step_size) + " | " + std::to_string(time) + " -> " + std::to_string(event_time));
     for (size_t i = 0; i < fmus.size(); ++i)
@@ -627,6 +666,8 @@ oms_status_enu_t oms::SystemSC::doStepEuler()
       if (oms_status_ok != status) return status;
     }
 
+    updateInputs(eventGraph);
+
     // b. Event Detection
     event_detected = event_time == tnext;
     logDebug("Event detected: " + std::to_string(event_detected));
@@ -657,10 +698,6 @@ oms_status_enu_t oms::SystemSC::doStepEuler()
         time = event_time;
         step_size_adjustment = maximumStepSize;
 
-        // set time
-        for (const auto& component : getComponents())
-          component.second->setTime(time);
-
         for (size_t i = 0; i < fmus.size(); ++i)
         {
           fmistatus = fmi2_completedIntegratorStep(fmus[i]->getFMU(), fmi2True, &callEventUpdate[i], &terminateSimulation[i]);
@@ -696,10 +733,6 @@ oms_status_enu_t oms::SystemSC::doStepEuler()
         if (isTopLevelSystem())
           getModel().emit(time, false);
 
-        // set time
-        for (const auto& component : getComponents())
-          component.second->setTime(time);
-
         // Enter event mode and handle discrete state updates for each FMU
         for (size_t i = 0; i < fmus.size(); ++i)
         {
@@ -713,7 +746,16 @@ oms_status_enu_t oms::SystemSC::doStepEuler()
 
           fmistatus = fmi2_enterContinuousTimeMode(fmus[i]->getFMU());
           if (fmi2OK != fmistatus) logError_FMUCall("fmi2_enterContinuousTimeMode", fmus[i]);
-
+        }
+        
+        updateInputs(eventGraph);
+        
+        // emit the right limit of the event
+        if (isTopLevelSystem())
+          getModel().emit(time, true);
+        
+        for (size_t i = 0; i < fmus.size(); ++i)
+        {
           if (nStates[i] > 0)
           {
             status = fmus[i]->getContinuousStates(states_backup[i]);
@@ -722,8 +764,11 @@ oms_status_enu_t oms::SystemSC::doStepEuler()
             if (oms_status_ok != status) return status;
           }
 
-          status = fmus[i]->getEventindicators(event_indicators_prev[i]);
-          if (oms_status_ok != status) return status;
+          if (nEventIndicators[i] > 0)
+          {
+            status = fmus[i]->getEventindicators(event_indicators_prev[i]);
+            if (oms_status_ok != status) return status;
+          }
         }
 
         // find next time event
@@ -740,11 +785,6 @@ oms_status_enu_t oms::SystemSC::doStepEuler()
             terminated = true;
           }
         }
-
-        // emit the right limit of the event
-        updateInputs(eventGraph);
-        if (isTopLevelSystem())
-          getModel().emit(time, true);
       }
       else
       {
@@ -766,13 +806,14 @@ oms_status_enu_t oms::SystemSC::doStepEuler()
   return oms_status_ok;
 }
 
-oms_status_enu_t oms::SystemSC::doStepCVODE()
+oms_status_enu_t oms::SystemSC::doStepCVODE(double stopTime)
 {
   fmi2Status fmistatus;
   oms_status_enu_t status;
-  int flag;
+  int flag = 0;
+  bool tnext_is_event = false;
 
-  const fmi2Real end_time = std::min(time + maximumStepSize, getModel().getStopTime());
+  const fmi2Real end_time = stopTime;
 
   // find next time event
   fmi2Real tnext = end_time+1.0;
@@ -788,19 +829,47 @@ oms_status_enu_t oms::SystemSC::doStepCVODE()
       time = end_time;
     }
   }
+
+  tnext_is_event = tnext <= end_time;
+  tnext = std::min(tnext, end_time);
+
   logDebug("tnext: " + std::to_string(tnext));
 
-  while (time < end_time)
-  {
+  //while (time < end_time)
+  //{
     logDebug("CVode: " + std::to_string(time) + " -> " + std::to_string(end_time));
     for (size_t j=0, k=0; j < fmus.size(); ++j)
       for (size_t i=0; i < nStates[j]; ++i, ++k)
         NV_Ith_S(solverData.cvode.y, k) = states[j][i];
 
-    flag = CVode(solverData.cvode.mem, std::min(tnext, end_time), solverData.cvode.y, &time, CV_NORMAL);
+    double cvode_time = tnext;
+    CVodeGetCurrentTime(solverData.cvode.mem, &cvode_time);
+    if (cvode_time < tnext) {
+      flag = CVode(solverData.cvode.mem, tnext, solverData.cvode.y, &cvode_time, CV_ONE_STEP);
+      if (flag < 0)
+        return logError("SUNDIALS_ERROR: CVode() failed with flag = " + std::to_string(flag));
+    }
+
+    if (cvode_time > tnext) {
+      // interpolate result to tnext (this shouldn't take any further steps)
+      flag = CVode(solverData.cvode.mem, tnext, solverData.cvode.y, &cvode_time, CV_NORMAL);
+      if (flag < 0)
+        return logError("SUNDIALS_ERROR: CVode() failed with flag = " + std::to_string(flag));
+    }
+
+    time = cvode_time;
+
+    CVodeGetErrWeights(solverData.cvode.mem, solverData.cvode.ewt);
+
+    // set time
+    for (const auto& component : getComponents())
+      component.second->setTime(time);
 
     for (size_t i = 0, j=0; i < fmus.size(); ++i)
     {
+      if (nStates[i] == 0)
+        continue;
+
       for (size_t k = 0; k < nStates[i]; k++, j++)
         states[i][k] = NV_Ith_S(solverData.cvode.y, j);
 
@@ -809,23 +878,25 @@ oms_status_enu_t oms::SystemSC::doStepCVODE()
       if (oms_status_ok != status) return status;
     }
 
-    if (flag == CV_ROOT_RETURN || time == tnext)
+    updateInputs(eventGraph);
+    if (isTopLevelSystem())
+      getModel().emit(time, false);
+    
+    for (size_t i = 0; i < fmus.size(); ++i)
     {
-      logDebug("event found!!! " + std::to_string(time));
-
-      // set time
-      for (const auto& component : getComponents())
-        component.second->setTime(time);
-
-      for (size_t i = 0; i < fmus.size(); ++i)
+      if (nStates[i] > 0)
       {
-        fmistatus = fmi2_completedIntegratorStep(fmus[i]->getFMU(), fmi2True, &callEventUpdate[i], &terminateSimulation[i]);
-        if (fmi2OK != fmistatus) return logError_FMUCall("fmi2_completedIntegratorStep", fmus[i]);
+        status = fmus[i]->getDerivatives(states_der[i]);
+        if (oms_status_ok != status) return status;
       }
+      
+      fmistatus = fmi2_completedIntegratorStep(fmus[i]->getFMU(), fmi2True, &callEventUpdate[i], &terminateSimulation[i]);
+      if (fmi2OK != fmistatus) return logError_FMUCall("fmi2_completedIntegratorStep", fmus[i]);
+    }
 
-      updateInputs(eventGraph);
-      if (isTopLevelSystem())
-        getModel().emit(time, false);
+    if (flag == CV_ROOT_RETURN || tnext_is_event && time == tnext)
+    {
+      logDebug("event found!!! " + std::to_string(time));
 
       // Enter event mode and handle discrete state updates for each FMU
       for (size_t i = 0; i < fmus.size(); ++i)
@@ -839,15 +910,6 @@ oms_status_enu_t oms::SystemSC::doStepCVODE()
         if (fmi2OK != fmistatus) logError_FMUCall("fmi2_enterContinuousTimeMode", fmus[i]);
       }
 
-      for (size_t i = 0; i < fmus.size(); ++i)
-      {
-        if (0 == nStates[i])
-          continue;
-
-        status = fmus[i]->getContinuousStates(states[i]);
-        if (oms_status_ok != status) return status;
-      }
-
       // find next time event
       tnext = end_time+1.0;
       for (size_t i = 0; i < fmus.size(); ++i)
@@ -862,6 +924,10 @@ oms_status_enu_t oms::SystemSC::doStepCVODE()
           time = end_time;
         }
       }
+
+      tnext_is_event = tnext <= end_time;
+      tnext = std::min(tnext, end_time);
+
       logDebug("tnext: " + std::to_string(tnext));
 
       // emit the right limit of the event
@@ -876,6 +942,8 @@ oms_status_enu_t oms::SystemSC::doStepCVODE()
 
         status = fmus[i]->getContinuousStates(states[i]);
         if (oms_status_ok != status) return status;
+        status = fmus[i]->getDerivatives(states_der[i]);
+        if (oms_status_ok != status) return status;
       }
 
       for (size_t j=0, k=0; j < fmus.size(); ++j)
@@ -885,36 +953,14 @@ oms_status_enu_t oms::SystemSC::doStepCVODE()
       flag = CVodeReInit(solverData.cvode.mem, time, solverData.cvode.y);
       if (flag < 0) return logError("SUNDIALS_ERROR: CVodeReInit() failed with flag = " + std::to_string(flag));
 
-      continue;
+      return oms_status_ok;
     }
 
-    if (flag == CV_SUCCESS)
-    {
+    if (flag >= 0) // Some kind of success
       logDebug("CVode completed successfully at t = " + std::to_string(time));
-
-      // set time
-      for (const auto& component : getComponents())
-        component.second->setTime(time);
-
-      for (size_t i = 0; i < fmus.size(); ++i)
-      {
-        fmistatus = fmi2_completedIntegratorStep(fmus[i]->getFMU(), fmi2True, &callEventUpdate[i], &terminateSimulation[i]);
-        if (fmi2OK != fmistatus) return logError_FMUCall("fmi2_completedIntegratorStep", fmus[i]);
-
-        if (0 == nStates[i])
-          continue;
-
-        status = fmus[i]->getContinuousStates(states[i]);
-        if (oms_status_ok != status) return status;
-      }
-
-      updateInputs(eventGraph);
-      if (isTopLevelSystem())
-        getModel().emit(time, false);
-    }
     else
       return logError("CVode failed with flag = " + std::to_string(flag));
-  }
+  //}
 
   return oms_status_ok;
 
@@ -930,14 +976,25 @@ oms_status_enu_t oms::SystemSC::stepUntil(double stopTime)
 
   // main simulation loop
   oms_status_enu_t status = oms_status_ok;
-  while (time < std::min(stopTime, getModel().getStopTime()) && oms_status_ok == status)
+  double endTime = std::min(stopTime, getModel().getStopTime());
+  double lastTime = time;
+  while (time < endTime && oms_status_ok == status)
   {
-    status = doStep();
+    if (solverMethod == oms_solver_sc_explicit_euler)
+      status = doStepEuler(endTime);
+    else if (solverMethod == oms_solver_sc_cvode)
+      status = doStepCVODE(endTime);
+    else
+      return logError_InternalError;
+
     if (status != oms_status_ok)
       logWarning("Bad return code at time " + std::to_string(time));
 
-    if (isTopLevelSystem() && Flags::ProgressBar())
+    if (isTopLevelSystem() && Flags::ProgressBar() && time > lastTime + maximumStepSize)
+    {
       Log::ProgressBar(startTime, stopTime, time);
+      lastTime = time;
+    }
   }
 
   if (isTopLevelSystem() && Flags::ProgressBar())
@@ -946,7 +1003,7 @@ oms_status_enu_t oms::SystemSC::stepUntil(double stopTime)
   return status;
 }
 
-oms_status_enu_t oms::SystemSC::updateInputs(DirectedGraph& graph)
+oms_status_enu_t oms::SystemSC::updateInputs(DirectedGraph& graph, double tolerance)
 {
   CallClock callClock(clock);
   oms_status_enu_t status;
@@ -994,7 +1051,7 @@ oms_status_enu_t oms::SystemSC::updateInputs(DirectedGraph& graph)
     }
     else
     {
-      status = solveAlgLoop(graph, loopNum);
+      status = solveAlgLoop(graph, loopNum, tolerance);
       if (oms_status_ok != status)
       {
         forceLoopsToBeUpdated();
diff --git a/src/OMSimulatorLib/SystemSC.h b/src/OMSimulatorLib/SystemSC.h
index f9b07d830..3a8cd3b2c 100644
--- a/src/OMSimulatorLib/SystemSC.h
+++ b/src/OMSimulatorLib/SystemSC.h
@@ -63,7 +63,7 @@ namespace oms
     oms_status_enu_t doStep();
     oms_status_enu_t stepUntil(double stopTime);
 
-    oms_status_enu_t updateInputs(DirectedGraph& graph);
+    oms_status_enu_t updateInputs(DirectedGraph& graph, double tolerance = 0.0);
 
     std::string getSolverName() const;
     oms_status_enu_t setSolverMethod(std::string);
@@ -71,8 +71,8 @@ namespace oms
     oms_status_enu_t setSolver(oms_solver_enu_t solver) {if (solver > oms_solver_sc_min && solver < oms_solver_sc_max) {solverMethod=solver; return oms_status_ok;} return oms_status_error;}
 
   private:
-    oms_status_enu_t doStepEuler();
-    oms_status_enu_t doStepCVODE();
+    oms_status_enu_t doStepEuler(double stopTime);
+    oms_status_enu_t doStepCVODE(double stopTime);
 
   protected:
     SystemSC(const ComRef& cref, Model* parentModel, System* parentSystem);
@@ -107,6 +107,7 @@ namespace oms
       SUNMatrix J;            /* Matrix used by linear solver */
       N_Vector liny;          /* Vector used by linear solver */
       N_Vector abstol;
+      N_Vector ewt;
     };
 
     union SolverData_t