diff --git a/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp b/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp
index deac97076c985c36f86ce8242296121fc82338f4..db9ef994f1c71e4a6cea64efb9f6fbff80251387 100644
--- a/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp
+++ b/Tests/NumLib/TestTimeSteppingEvolutionaryPIDcontroller.cpp
@@ -28,7 +28,7 @@ std::unique_ptr<NumLib::TimeStepAlgorithm> createTestTimeStepper(
     BaseLib::ConfigTree conf(std::move(ptree), "", BaseLib::ConfigTree::onerror,
                              BaseLib::ConfigTree::onwarning);
     auto const& sub_config = conf.getConfigSubtree("time_stepping");
-    return NumLib::createEvolutionaryPIDcontroller(sub_config);
+    return NumLib::createEvolutionaryPIDcontroller(sub_config, {});
 }
 
 TEST(NumLibTimeStepping, testEvolutionaryPIDcontroller)
diff --git a/Tests/NumLib/TestTimeSteppingFixed.cpp b/Tests/NumLib/TestTimeSteppingFixed.cpp
index 18d67fc806af3652e2df19e594ffc1c6e761a6d5..a982f0c21a82247a312bf6849c64355cb444c8c0 100644
--- a/Tests/NumLib/TestTimeSteppingFixed.cpp
+++ b/Tests/NumLib/TestTimeSteppingFixed.cpp
@@ -12,13 +12,22 @@
 
 #include <gtest/gtest.h>
 
+#include <numeric>
 #include <vector>
 
+#include "NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h"
 #include "NumLib/TimeStepping/Algorithms/FixedTimeStepping.h"
 #include "NumLib/TimeStepping/TimeStep.h"
 #include "Tests/TestTools.h"
 #include "TimeSteppingTestingTools.h"
 
+namespace NumLib
+{
+extern void incorporateFixedTimesForOutput(
+    double t_initial, double t_end, std::vector<double>& timesteps,
+    std::vector<double> const& fixed_times_for_output);
+}
+
 TEST(NumLib, TimeSteppingFixed)
 {
     std::vector<int> const dummy_number_iterations = {0, 0, 0, 0};
@@ -77,3 +86,90 @@ TEST(NumLib, TimeSteppingFixed)
                           std::numeric_limits<double>::epsilon());
     }
 }
+
+TEST(NumLib, TimeSteppingFixed_incorporateFixedTimesForOutput)
+{
+    double t_initial = 1.0;
+    std::vector<double> timesteps{{10, 10, 10}};
+    std::vector<double> fixed_times_for_output{{9, 12, 28}};
+
+    auto const expected_time =
+        std::accumulate(timesteps.begin(), timesteps.end(), 0.0);
+
+    NumLib::incorporateFixedTimesForOutput(t_initial, expected_time, dts,
+                                           fixed_times_for_output);
+
+    ASSERT_EQ(expected_time,
+              std::accumulate(timesteps.begin(), timesteps.end(), 0.0));
+    ASSERT_EQ(8.0, timesteps[0]);
+    ASSERT_EQ(2.0, timesteps[1]);
+    ASSERT_EQ(1.0, timesteps[2]);
+    ASSERT_EQ(9.0, timesteps[3]);
+    ASSERT_EQ(7.0, timesteps[4]);
+    ASSERT_EQ(3.0, timesteps[5]);
+}
+
+TEST(NumLib, TimeSteppingFixed_incorporateFixedTimesForOutput_Matching)
+{
+    double t_initial = 1.0;
+    std::vector<double> timesteps{{10, 10, 10}};
+    std::vector<double> fixed_times_for_output{{9, 11, 31}};
+
+    auto const expected_time =
+        std::accumulate(timesteps.begin(), timesteps.end(), 0.0);
+
+    NumLib::incorporateFixedTimesForOutput(t_initial, expected_time, timesteps,
+                                           fixed_times_for_output);
+
+    ASSERT_EQ(expected_time,
+              std::accumulate(timesteps.begin(), timesteps.end(), 0.0));
+    ASSERT_EQ(8.0, timesteps[0]);
+    ASSERT_EQ(2.0, timesteps[1]);
+    ASSERT_EQ(10.0, timesteps[2]);
+    ASSERT_EQ(10.0, timesteps[3]);
+}
+
+TEST(
+    NumLib,
+    TimeSteppingFixed_incorporateFixedTimesForOutput_OutputTimeBeforeSimulationStartTime)
+{
+    double t_initial = 10.0;
+    std::vector<double> timesteps{{10, 10, 10}};
+    std::vector<double> fixed_times_for_output{{9, 12, 28}};
+
+    auto const expected_time =
+        std::accumulate(timesteps.begin(), timesteps.end(), t_initial);
+
+    NumLib::incorporateFixedTimesForOutput(t_initial, expected_time, timesteps,
+                                           fixed_times_for_output);
+    ASSERT_EQ(expected_time,
+              std::accumulate(timesteps.begin(), timesteps.end(), t_initial));
+    ASSERT_EQ(2.0, timesteps[0]);
+    ASSERT_EQ(8.0, timesteps[1]);
+    ASSERT_EQ(8.0, timesteps[2]);
+    ASSERT_EQ(2.0, timesteps[3]);
+    ASSERT_EQ(10.0, timesteps[4]);
+}
+
+TEST(
+    NumLib,
+    TimeSteppingFixed_incorporateFixedTimesForOutput_OutputTimeAfterSimulationEndTime)
+{
+    double t_initial = 1.0;
+    std::vector<double> timesteps{{10, 10, 10}};
+    std::vector<double> fixed_times_for_output{{9, 12, 28, 33}};
+
+    auto const expected_time =
+        std::accumulate(timesteps.begin(), timesteps.end(), t_initial);
+
+    NumLib::incorporateFixedTimesForOutput(t_initial, expected_time, timesteps,
+                                           fixed_times_for_output);
+    ASSERT_EQ(expected_time,
+              std::accumulate(timesteps.begin(), timesteps.end(), t_initial));
+    ASSERT_EQ(8.0, timesteps[0]);
+    ASSERT_EQ(2.0, timesteps[1]);
+    ASSERT_EQ(1.0, timesteps[2]);
+    ASSERT_EQ(9.0, timesteps[3]);
+    ASSERT_EQ(7.0, timesteps[4]);
+    ASSERT_EQ(3.0, timesteps[5]);
+}
diff --git a/Tests/NumLib/TestTimeSteppingIterationNumber.cpp b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
index e7681be237abb30efe9bc0cebb35c51af64f4e87..5de4511128140f2e0e071170bdc560fc9b37dbb1 100644
--- a/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
+++ b/Tests/NumLib/TestTimeSteppingIterationNumber.cpp
@@ -25,9 +25,9 @@ TEST(NumLib, TimeSteppingIterationNumberBased1)
 {
     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,
-                                                 std::move(iter_times_vector),
-                                                 std::move(multiplier_vector));
+    NumLib::IterationNumberBasedTimeStepping alg(
+        1, 31, 1, 10, 1, std::move(iter_times_vector),
+        std::move(multiplier_vector), {});
 
     const double solution_error = 0.;
     const double end_time = alg.end();
@@ -145,9 +145,9 @@ TEST(NumLib, TimeSteppingIterationNumberBased2)
 {
     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,
-                                                 std::move(iter_times_vector),
-                                                 std::move(multiplier_vector));
+    NumLib::IterationNumberBasedTimeStepping alg(
+        1, 31, 1, 10, 1, std::move(iter_times_vector),
+        std::move(multiplier_vector), {});
 
     std::vector<int> nr_iterations = {0, 2, 2, 2, 4, 6, 8, 4, 1};
     const std::vector<double> expected_vec_t = {1,  2,  4,  8,  16,
@@ -165,9 +165,9 @@ TEST(NumLib, TimeSteppingIterationNumberBased2FixedOutputTimes)
     std::vector<int> iter_times_vector = {0, 3, 5, 7};
     std::vector<double> multiplier_vector = {2.0, 1.0, 0.5, 0.25};
     std::vector<double> fixed_output_times = {5, 20};
-    NumLib::IterationNumberBasedTimeStepping alg(1, 31, 1, 10, 1,
-                                                 std::move(iter_times_vector),
-                                                 std::move(multiplier_vector));
+    NumLib::IterationNumberBasedTimeStepping alg(
+        1, 31, 1, 10, 1, std::move(iter_times_vector),
+        std::move(multiplier_vector), {});
 
     std::vector<int> nr_iterations = {0, 2, 2, 2, 4, 6, 8, 4, 1, 1, 1, 1, 1};
     const std::vector<double> expected_vec_t = {1,  2,  4,  5,  7,  9,  10,