diff --git a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp
index ff891e6503b96805515dbeb94f5a5d6884f703fd..985fb0a4f89763264561ccc2936ca63ef361cbed 100644
--- a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp
+++ b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.cpp
@@ -108,9 +108,10 @@ setInitialConditions(ProjectData& project,
 bool
 UncoupledProcessesTimeLoop::
 solveOneTimeStepOneProcess(
-        GlobalVector& x, double const t, double const delta_t,
+        GlobalVector& x, std::size_t timestep, double const t, double const delta_t,
         SingleProcessData& process_data,
-        UncoupledProcessesTimeLoop::Process& process)
+        UncoupledProcessesTimeLoop::Process& process,
+        ProcessLib::Output const& output_control)
 {
     auto& time_disc        =  process.getTimeDiscretization();
     auto& ode_sys          = *process_data.tdisc_ode_sys;
@@ -130,7 +131,13 @@ solveOneTimeStepOneProcess(
 
     process.preTimestep(x, t, delta_t);
 
-    bool nonlinear_solver_succeeded = nonlinear_solver.solve(x);
+    auto const post_iteration_callback = [&](unsigned iteration,
+                                             GlobalVector const& x) {
+        output_control.doOutputIteration(process, timestep, t, x, iteration);
+    };
+
+    bool nonlinear_solver_succeeded =
+        nonlinear_solver.solve(x, post_iteration_callback);
 
     auto& mat_strg = process_data.mat_strg;
     time_disc.pushState(t, x, mat_strg);
@@ -186,7 +193,8 @@ bool UncoupledProcessesTimeLoop::loop(ProjectData& project)
             auto& x = *_process_solutions[pcs_idx];
 
             nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
-                        x, t, delta_t, per_process_data[pcs_idx], **p);
+                x, timestep, t, delta_t, per_process_data[pcs_idx], **p,
+                out_ctrl);
 
             if (!nonlinear_solver_succeeded) {
                 ERR("The nonlinear solver failed in timestep #%u at t = %g s"
diff --git a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
index a89d8d15826e6f1b28df181a68ef95cb01387660..2a9435f2eb4b81ae2b5576a59ad85930682c2836 100644
--- a/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
+++ b/Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h
@@ -155,9 +155,9 @@ private:
 
     //! Solves one timestep for the given \c process.
     bool solveOneTimeStepOneProcess(
-            GlobalVector& x, double const t, double const delta_t,
+            GlobalVector& x, std::size_t timestep, double const t, double const delta_t,
             SingleProcessData& process_data,
-            Process& process);
+            Process& process, ProcessLib::Output const& output_control);
 
     //! Sets the EquationSystem for the given nonlinear solver,
     //! which is Picard or Newton depending on the NLTag.
diff --git a/NumLib/ODESolver/NonlinearSolver-impl.h b/NumLib/ODESolver/NonlinearSolver-impl.h
index 418bcea51048374e2c23a806c69ffa84827b8c93..8b1fb2ece36056a9740ff3eab5b9c5086bdf81df 100644
--- a/NumLib/ODESolver/NonlinearSolver-impl.h
+++ b/NumLib/ODESolver/NonlinearSolver-impl.h
@@ -32,10 +32,10 @@ assemble(Vector const& x) const
     _equation_system->assembleMatricesPicard(x);
 }
 
-template<typename Matrix, typename Vector>
-bool
-NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard>::
-solve(Vector &x)
+template <typename Matrix, typename Vector>
+bool NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Picard>::solve(
+    Vector& x,
+    std::function<void(unsigned, Vector const&)> const& postIterationCallback)
 {
     namespace BLAS = MathLib::BLAS;
     auto& sys = *_equation_system;
@@ -71,6 +71,9 @@ solve(Vector &x)
         }
         else
         {
+            if (postIterationCallback)
+                postIterationCallback(iteration, x_new);
+
             switch(sys.postIteration(x_new))
             {
             case IterationResult::SUCCESS:
@@ -143,10 +146,10 @@ assemble(Vector const& x) const
     //      equation every time and could not forget it.
 }
 
-template<typename Matrix, typename Vector>
-bool
-NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::
-solve(Vector &x)
+template <typename Matrix, typename Vector>
+bool NonlinearSolver<Matrix, Vector, NonlinearSolverTag::Newton>::solve(
+    Vector& x,
+    std::function<void(unsigned, Vector const&)> const& postIterationCallback)
 {
     namespace BLAS = MathLib::BLAS;
     auto& sys = *_equation_system;
@@ -195,6 +198,9 @@ solve(Vector &x)
                     MathLib::GlobalVectorProvider<Vector>::provider.getVector(x, _x_new_id);
             BLAS::axpy(x_new, -_alpha, minus_delta_x);
 
+            if (postIterationCallback)
+                postIterationCallback(iteration, x_new);
+
             switch(sys.postIteration(x_new))
             {
             case IterationResult::SUCCESS:
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index aa4d4cad51170018c4790694ab171e30d370dbfb..aaed26ca2778b84b875c0c354a8041061df49ec2 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -53,11 +53,14 @@ public:
     /*! Assemble and solve the equation system.
      *
      * \param x   in: the initial guess, out: the solution.
+     * \param postIterationCallback called after each iteration if set.
      *
      * \retval true if the equation system could be solved
      * \retval false otherwise
      */
-    virtual bool solve(Vector& x) = 0;
+    virtual bool solve(Vector& x,
+                       std::function<void(unsigned, Vector const&)> const&
+                           postIterationCallback) = 0;
 
     virtual ~NonlinearSolverBase() = default;
 };
@@ -108,7 +111,9 @@ public:
 
     void assemble(Vector const& x) const override;
 
-    bool solve(Vector& x) override;
+    bool solve(Vector& x,
+               std::function<void(unsigned, Vector const&)> const&
+                   postIterationCallback) override;
 
 private:
     LinearSolver& _linear_solver;
@@ -159,7 +164,9 @@ public:
 
     void assemble(Vector const& x) const override;
 
-    bool solve(Vector& x) override;
+    bool solve(Vector& x,
+               std::function<void(unsigned, Vector const&)> const&
+                   postIterationCallback) override;
 
 private:
     LinearSolver& _linear_solver;
diff --git a/NumLib/ODESolver/TimeLoopSingleODE.h b/NumLib/ODESolver/TimeLoopSingleODE.h
index f7c45ad92ce0e61b2263cebbb11efca5ad1690fb..a974e1ba12a9883f80344fc387d91fde9e8eb1df 100644
--- a/NumLib/ODESolver/TimeLoopSingleODE.h
+++ b/NumLib/ODESolver/TimeLoopSingleODE.h
@@ -110,7 +110,7 @@ loop(const double t0, Vector const& x0, const double t_end, const double delta_t
         // INFO("time: %e, delta_t: %e", t, delta_t);
         time_disc.nextTimestep(t, delta_t);
 
-        nl_slv_succeeded = _nonlinear_solver->solve(x);
+        nl_slv_succeeded = _nonlinear_solver->solve(x, nullptr);
         if (!nl_slv_succeeded) break;
 
         time_disc.pushState(t, x, _ode_sys);