diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp
index 17a1a7dfc392a5b3118ee183358eca1635759ef8..ef6566cc0a384a288cf499cd3820cf054c334cf7 100644
--- a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp
@@ -62,7 +62,6 @@ bool IterationNumberBasedTimeStepping::next(double const /*solution_error*/,
     // prepare the next time step info
     _ts_current = _ts_prev;
     _ts_current += getNextTimeStepSize();
-    _iter_times = 0;
 
     return true;
 }
diff --git a/Tests/NumLib/TestODEInt.cpp b/Tests/NumLib/TestODEInt.cpp
index aa110efc8e4afa2eb0abb31cbc77b8e13f19b29d..bc79134f9bbcbaba3ad222ad8539770f9db95777 100644
--- a/Tests/NumLib/TestODEInt.cpp
+++ b/Tests/NumLib/TestODEInt.cpp
@@ -104,7 +104,7 @@ public:
 
         if (num_timesteps > 0)
         {
-            EXPECT_TRUE(loop.loop(t0, x0, t_end, delta_t, cb));
+            EXPECT_TRUE(loop.loop(t0, x0, t_end, delta_t, cb).error_norms_met);
         }
 
         for (auto& x :  sol.solutions)
diff --git a/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp b/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp
index 45b1c69d0288e7cb86f55ed92c7b77e98065d26b..debb22fd2203d30bd7e005390768eda2ffbf9dae 100644
--- a/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp
+++ b/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp
@@ -50,8 +50,9 @@ TEST(NumLibTimeStepping, testEvolutionaryPIDcontroller)
     auto const PIDStepper = createTestTimeStepper(xml);
 
     double solution_error = 0.;
+    int const number_iterations = 0;
     // 1st step
-    ASSERT_TRUE(PIDStepper->next(solution_error));
+    ASSERT_TRUE(PIDStepper->next(solution_error, number_iterations));
     NumLib::TimeStep ts = PIDStepper->getTimeStep();
     double h_new = 0.01;
     double t_previous = 0.;
@@ -64,7 +65,7 @@ TEST(NumLibTimeStepping, testEvolutionaryPIDcontroller)
 
     // e_n_minus1 is filled.
     solution_error = 1.0e-4;
-    PIDStepper->next(solution_error);
+    PIDStepper->next(solution_error, number_iterations);
     ts = PIDStepper->getTimeStep();
     h_new = ts.dt();
     ASSERT_EQ(2u, ts.steps());
@@ -76,7 +77,7 @@ TEST(NumLibTimeStepping, testEvolutionaryPIDcontroller)
 
     // e_n_minus2 is filled.
     solution_error = 0.5e-3;
-    PIDStepper->next(solution_error);
+    PIDStepper->next(solution_error, number_iterations);
     ts = PIDStepper->getTimeStep();
     h_new = ts.dt();
     ASSERT_EQ(3u, ts.steps());
@@ -86,7 +87,7 @@ TEST(NumLibTimeStepping, testEvolutionaryPIDcontroller)
 
     // error > TOL=1.3-3, step rejected and new step size estimated.
     solution_error = 0.01;
-    PIDStepper->next(solution_error);
+    PIDStepper->next(solution_error, number_iterations);
     ts = PIDStepper->getTimeStep();
     h_new = ts.dt();
     // No change in ts.steps
@@ -100,7 +101,7 @@ TEST(NumLibTimeStepping, testEvolutionaryPIDcontroller)
 
     // With e_n, e_n_minus1, e_n_minus2
     solution_error = 0.4e-3;
-    PIDStepper->next(solution_error);
+    PIDStepper->next(solution_error, number_iterations);
     ts = PIDStepper->getTimeStep();
     h_new = ts.dt();
     ASSERT_EQ(4u, ts.steps());
diff --git a/Tests/NumLib/TestTimeSteppingFixed.cpp b/Tests/NumLib/TestTimeSteppingFixed.cpp
index 7e462f7e38555770fc86cff36f2a67c47bcd3e2f..c1caeaec8f41c10a76d2a779c7134f78ac43072e 100644
--- a/Tests/NumLib/TestTimeSteppingFixed.cpp
+++ b/Tests/NumLib/TestTimeSteppingFixed.cpp
@@ -15,23 +15,26 @@
 
 #include <logog/include/logog.hpp>
 
-#include "NumLib/TimeStepping/TimeStep.h"
 #include "NumLib/TimeStepping/Algorithms/FixedTimeStepping.h"
+#include "NumLib/TimeStepping/TimeStep.h"
 
 #include "Tests/TestTools.h"
 #include "TimeSteppingTestingTools.h"
 
 TEST(NumLib, TimeSteppingFixed)
 {
+    std::vector<int> const dummy_number_iterations = {0, 0, 0, 0};
     // homogeneous dt
     {
         NumLib::FixedTimeStepping fixed(1, 31, 10);
         const std::vector<double> expected_vec_t = {1, 11, 21, 31};
 
-        std::vector<double> vec_t = timeStepping(fixed);
+        std::vector<double> vec_t =
+            timeStepping(fixed, dummy_number_iterations);
 
         ASSERT_EQ(expected_vec_t.size(), vec_t.size());
-        ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), std::numeric_limits<double>::epsilon());
+        ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(),
+                          std::numeric_limits<double>::epsilon());
     }
 
     // dt vector (t_end == t0 + sum(dt))
@@ -40,10 +43,12 @@ TEST(NumLib, TimeSteppingFixed)
         NumLib::FixedTimeStepping fixed(1, 31, 10);
         const std::vector<double> expected_vec_t = {1, 11, 21, 31};
 
-        std::vector<double> vec_t = timeStepping(fixed);
+        std::vector<double> vec_t =
+            timeStepping(fixed, dummy_number_iterations);
 
         ASSERT_EQ(expected_vec_t.size(), vec_t.size());
-        ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), std::numeric_limits<double>::epsilon());
+        ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(),
+                          std::numeric_limits<double>::epsilon());
     }
 
     // dt vector (t_end < t0 + sum(dt))
@@ -52,10 +57,12 @@ TEST(NumLib, TimeSteppingFixed)
         NumLib::FixedTimeStepping fixed(1, 31, fixed_dt);
         const std::vector<double> expected_vec_t = {1, 6, 16, 31};
 
-        std::vector<double> vec_t = timeStepping(fixed);
+        std::vector<double> vec_t =
+            timeStepping(fixed, dummy_number_iterations);
 
         ASSERT_EQ(expected_vec_t.size(), vec_t.size());
-        ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), std::numeric_limits<double>::epsilon());
+        ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(),
+                          std::numeric_limits<double>::epsilon());
     }
 
     // dt vector (t_end > t0 + sum(dt))
@@ -64,9 +71,11 @@ TEST(NumLib, TimeSteppingFixed)
         NumLib::FixedTimeStepping fixed(1, 31, fixed_dt);
         const std::vector<double> expected_vec_t = {1, 6, 16, 26};
 
-        std::vector<double> vec_t = timeStepping(fixed);
+        std::vector<double> vec_t =
+            timeStepping(fixed, dummy_number_iterations);
 
         ASSERT_EQ(expected_vec_t.size(), vec_t.size());
-        ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(), std::numeric_limits<double>::epsilon());
+        ASSERT_ARRAY_NEAR(expected_vec_t, vec_t, expected_vec_t.size(),
+                          std::numeric_limits<double>::epsilon());
     }
 }
diff --git a/Tests/NumLib/TestTimeSteppingIterationNumber.cpp b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
index 967588793920e24d06c4518b833180f070c61cd8..0f29a6ad6a80b4e2b257559922fd99e1e4743ace 100644
--- a/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
+++ b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
@@ -24,13 +24,15 @@
 
 TEST(NumLib, TimeSteppingIterationNumberBased1)
 {
-    std::vector<std::size_t> iter_times_vector = {0, 3, 5, 7};
+    std::vector<int> iter_times_vector = {0, 3, 5, 7};
     std::vector<double> multiplier_vector = {2.0, 1.0, 0.5, 0.25};
-    NumLib::IterationNumberBasedTimeStepping alg(1, 31, 1, 10, 1, iter_times_vector, multiplier_vector);
+    NumLib::IterationNumberBasedTimeStepping alg(1, 31, 1, 10, 1,
+                                                 std::move(iter_times_vector),
+                                                 std::move(multiplier_vector));
 
     const double solution_error = 0.;
 
-    ASSERT_TRUE(alg.next(solution_error));  // t=2, dt=1
+    ASSERT_TRUE(alg.next(solution_error, 1));  // t=2, dt=1
     NumLib::TimeStep ts = alg.getTimeStep();
     ASSERT_EQ(1u, ts.steps());
     ASSERT_EQ(1., ts.previous());
@@ -38,11 +40,10 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
     ASSERT_EQ(1., ts.dt());
     ASSERT_TRUE(alg.accepted());
 
-    ASSERT_TRUE(alg.next(solution_error));  // t=4, dt=2
+    ASSERT_TRUE(alg.next(solution_error, 1));  // t=4, dt=2
 
     // dt*=2
-    alg.setIterationNumber(3);
-    ASSERT_TRUE(alg.next(solution_error));  // t=8, dt=4
+    ASSERT_TRUE(alg.next(solution_error, 3));  // t=8, dt=4
     ts = alg.getTimeStep();
     ASSERT_EQ(3u, ts.steps());
     ASSERT_EQ(4., ts.previous());
@@ -51,8 +52,7 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
     ASSERT_TRUE(alg.accepted());
 
     // dt*=1
-    alg.setIterationNumber(5);
-    ASSERT_TRUE(alg.next(solution_error));  // t=12, dt=4
+    ASSERT_TRUE(alg.next(solution_error, 5));  // t=12, dt=4
     ts = alg.getTimeStep();
     ASSERT_EQ(4u, ts.steps());
     ASSERT_EQ(8., ts.previous());
@@ -61,8 +61,7 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
     ASSERT_TRUE(alg.accepted());
 
     // dt*=0.5
-    alg.setIterationNumber(7);
-    ASSERT_TRUE(alg.next(solution_error));  // t=14, dt=2
+    ASSERT_TRUE(alg.next(solution_error, 7));  // t=14, dt=2
     ts = alg.getTimeStep();
     ASSERT_EQ(5u, ts.steps());
     ASSERT_EQ(12., ts.previous());
@@ -71,8 +70,8 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
     ASSERT_TRUE(alg.accepted());
 
     // dt*=0.25 but dt_min = 1
-    alg.setIterationNumber(8);              // exceed max
-    ASSERT_TRUE(alg.next(solution_error));  // t=13, dt=1
+    ASSERT_TRUE(
+        alg.next(solution_error, 8 /* exceed maximum */));  // t=13, dt=1
     ts = alg.getTimeStep();
     ASSERT_EQ(5u, ts.steps());
     ASSERT_EQ(12., ts.previous());
@@ -81,8 +80,7 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
     ASSERT_FALSE(alg.accepted());
 
     // restart, dt*=1
-    alg.setIterationNumber(4);
-    ASSERT_TRUE(alg.next(solution_error));  // t=14, dt=1
+    ASSERT_TRUE(alg.next(solution_error, 4));  // t=14, dt=1
     ts = alg.getTimeStep();
     ASSERT_EQ(6u, ts.steps());
     ASSERT_EQ(13., ts.previous());
@@ -93,22 +91,23 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
 
 TEST(NumLib, TimeSteppingIterationNumberBased2)
 {
-    std::vector<std::size_t> iter_times_vector = {0, 3, 5, 7};
+    std::vector<int> iter_times_vector = {0, 3, 5, 7};
     std::vector<double> multiplier_vector = {2.0, 1.0, 0.5, 0.25};
-    NumLib::IterationNumberBasedTimeStepping alg(1, 31, 1, 10, 1, iter_times_vector, multiplier_vector);
+    NumLib::IterationNumberBasedTimeStepping alg(1, 31, 1, 10, 1,
+                                                 std::move(iter_times_vector),
+                                                 std::move(multiplier_vector));
 
-    std::vector<std::size_t> nr_iterations = {2, 2, 2, 4, 6, 8, 4, 4, 2, 2};
+    std::vector<int> nr_iterations = {2, 2, 2, 4, 6, 8, 4, 4, 2, 2};
     const std::vector<double> expected_vec_t = {1, 2, 4, 8, 16, 24, 26, 28, 30, 31};
 
     struct IterationNumberUpdate
     {
-        IterationNumberUpdate(std::vector<std::size_t> vec,
-                              std::size_t& counter)
+        IterationNumberUpdate(std::vector<int> vec, std::size_t& counter)
             : _nr_iterations(std::move(vec)), i(counter)
         {
         }
 
-        std::vector<std::size_t> _nr_iterations;
+        std::vector<int> _nr_iterations;
         std::size_t& i;
 
         void operator()(NumLib::IterationNumberBasedTimeStepping &obj)
@@ -121,7 +120,7 @@ TEST(NumLib, TimeSteppingIterationNumberBased2)
     std::size_t counter = 0;
     IterationNumberUpdate update(nr_iterations, counter);
 
-    std::vector<double> vec_t = timeStepping(alg, &update);
+    std::vector<double> vec_t = timeStepping(alg, nr_iterations, &update);
 
     ASSERT_EQ(expected_vec_t.size(), vec_t.size());
     ASSERT_EQ(1u, alg.getNumberOfRepeatedSteps());
diff --git a/Tests/NumLib/TimeLoopSingleODE.h b/Tests/NumLib/TimeLoopSingleODE.h
index 31dcf3809abddf8cecd44159f800b2331f3fb03d..5b370a6feb50faba7c7ee56c7880ef17d6a88dbf 100644
--- a/Tests/NumLib/TimeLoopSingleODE.h
+++ b/Tests/NumLib/TimeLoopSingleODE.h
@@ -66,8 +66,9 @@ public:
      * \retval false otherwise
      */
     template <typename Callback>
-    bool loop(const double t0, GlobalVector const& x0, const double t_end,
-              const double delta_t, Callback& post_timestep);
+    NumLib::NonlinearSolverStatus loop(const double t0, GlobalVector const& x0,
+                                       const double t_end, const double delta_t,
+                                       Callback& post_timestep);
 
 private:
     TDiscODESys& _ode_sys;
@@ -80,13 +81,12 @@ private:
 
 template <NonlinearSolverTag NLTag>
 template <typename Callback>
-bool TimeLoopSingleODE<NLTag>::loop(const double t0, GlobalVector const& x0,
-                                    const double t_end, const double delta_t,
-                                    Callback& post_timestep)
+NumLib::NonlinearSolverStatus TimeLoopSingleODE<NLTag>::loop(
+    const double t0, GlobalVector const& x0, const double t_end,
+    const double delta_t, Callback& post_timestep)
 {
     // solution vector
-    GlobalVector& x =
-        NumLib::GlobalVectorProvider::provider.getVector(x0);
+    GlobalVector& x = NumLib::GlobalVectorProvider::provider.getVector(x0);
 
     auto& time_disc = _ode_sys.getTimeDiscretization();
 
@@ -104,7 +104,7 @@ bool TimeLoopSingleODE<NLTag>::loop(const double t0, GlobalVector const& x0,
 
     double t;
     unsigned timestep = 0;
-    bool nl_slv_succeeded = true;
+    NumLib::NonlinearSolverStatus nonlinear_solver_status = {false, 0};
     for (t = t0 + delta_t; t < t_end + std::numeric_limits<double>::epsilon();
          t = t0 + (timestep + 1) * delta_t)
     {
@@ -113,8 +113,8 @@ bool TimeLoopSingleODE<NLTag>::loop(const double t0, GlobalVector const& x0,
         // INFO("time: %e, delta_t: %e", t, delta_t);
         time_disc.nextTimestep(t, delta_t);
 
-        nl_slv_succeeded = _nonlinear_solver->solve(x, nullptr);
-        if (!nl_slv_succeeded)
+        nonlinear_solver_status = _nonlinear_solver->solve(x, nullptr);
+        if (!nonlinear_solver_status.error_norms_met)
             break;
 
         time_disc.pushState(t, x, _ode_sys);
@@ -127,10 +127,10 @@ bool TimeLoopSingleODE<NLTag>::loop(const double t0, GlobalVector const& x0,
 
     NumLib::GlobalVectorProvider::provider.releaseVector(x);
 
-    if (!nl_slv_succeeded)
+    if (!nonlinear_solver_status.error_norms_met)
     {
         ERR("Nonlinear solver failed in timestep #%u at t = %g s", timestep, t);
     }
-    return nl_slv_succeeded;
-}
+    return nonlinear_solver_status;
 }
+}  // namespace NumLib
diff --git a/Tests/NumLib/TimeSteppingTestingTools.h b/Tests/NumLib/TimeSteppingTestingTools.h
index 00e8986445489dcbc8da29c82f870facd46f329d..b499a760d421618391ea6c68b107f0498533bb47 100644
--- a/Tests/NumLib/TimeSteppingTestingTools.h
+++ b/Tests/NumLib/TimeSteppingTestingTools.h
@@ -25,14 +25,22 @@ struct Dummy
     void operator()(T &/*obj*/) {}
 };
 
-template <class T_TIME_STEPPING, class T=Dummy>
-std::vector<double> timeStepping(T_TIME_STEPPING &algorithm, T* obj=nullptr)
+template <class T_TIME_STEPPING, class T = Dummy>
+std::vector<double> timeStepping(T_TIME_STEPPING& algorithm,
+                                 std::vector<int> const& number_iterations,
+                                 T* obj = nullptr)
 {
     std::vector<double> vec_t;
     vec_t.push_back(algorithm.begin());
 
-    while (algorithm.next(0.0))
+    double const solution_error = 0;
+    for (auto const& i : number_iterations)
     {
+        if (!algorithm.next(solution_error, i))
+        {
+            break;
+        }
+
         NumLib::TimeStep t = algorithm.getTimeStep();
         //INFO("t: n=%d,t=%g,dt=%g", t.steps(), t.current(), t.dt());
         if (obj)