diff --git a/src/solve/Solver.cpp b/src/solve/Solver.cpp
index 6ac9ea21e..7f9bf06a7 100644
--- a/src/solve/Solver.cpp
+++ b/src/solve/Solver.cpp
@@ -56,8 +56,8 @@ Solver::Solver(const nlohmann::json& config) {
   sanity_checks();
 }
 
-void Solver::run() {
-  auto state = initial_state;
+void Solver::setup_initial() {
+  state = initial_state;
 
   // Create steady initial condition
   if (simparams.sim_steady_initial) {
@@ -79,12 +79,15 @@ void Solver::run() {
   // Use the initial condition (steady or user-provided) to set up parameters
   // which depend on the initial condition
   this->model->setup_initial_state_dependent_parameters(state);
+  DEBUG_MSG("Setup initial");
+}
 
+void Solver::setup_integrator() {
   // Set-up integrator
   DEBUG_MSG("Setup time integration");
-  Integrator integrator(this->model.get(), simparams.sim_time_step_size,
-                        simparams.sim_rho_infty, simparams.sim_abs_tol,
-                        simparams.sim_nliter);
+  integrator = Integrator(this->model.get(), simparams.sim_time_step_size,
+                          simparams.sim_rho_infty, simparams.sim_abs_tol,
+                          simparams.sim_nliter);
 
   // Initialize loop
   states = std::vector<State>();
@@ -102,8 +105,10 @@ void Solver::run() {
     states.reserve(num_states);
     times.reserve(num_states);
   }
-  double time = 0.0;
+  time = 0.0;
+}
 
+void Solver::run_integration() {
   // Run integrator
   DEBUG_MSG("Run time integration");
   int interval_counter = 0;
@@ -113,31 +118,37 @@ void Solver::run() {
   if (simparams.output_all_cycles || (0 >= start_last_cycle)) {
     times.push_back(time);
     states.push_back(std::move(state));
+    DEBUG_MSG("Added initial state and time");
   }
 
   int num_time_pts_in_two_cycles;
   std::vector<State> states_last_two_cycles;
   int last_two_cycles_time_pt_counter = 0;
+
   if (simparams.use_cycle_to_cycle_error) {
     num_time_pts_in_two_cycles = 2 * (simparams.sim_pts_per_cycle - 1) + 1;
     states_last_two_cycles =
         std::vector<State>(num_time_pts_in_two_cycles, state);
+    DEBUG_MSG("Initialized cycle to cycle error tracking with "
+              << num_time_pts_in_two_cycles << " points");
   }
+
   for (int i = 1; i < simparams.sim_num_time_steps; i++) {
     if (simparams.use_cycle_to_cycle_error) {
       if (i == simparams.sim_num_time_steps - num_time_pts_in_two_cycles + 1) {
-        // add first state
         states_last_two_cycles[last_two_cycles_time_pt_counter] = state;
-        last_two_cycles_time_pt_counter +=
-            1;  // last_two_cycles_time_pt_counter becomes 1
+        last_two_cycles_time_pt_counter += 1;
       }
     }
+
     state = integrator.step(state, time);
+
     if (simparams.use_cycle_to_cycle_error &&
         last_two_cycles_time_pt_counter > 0) {
       states_last_two_cycles[last_two_cycles_time_pt_counter] = state;
       last_two_cycles_time_pt_counter += 1;
     }
+
     interval_counter += 1;
     time = simparams.sim_time_step_size * double(i);
 
@@ -166,9 +177,11 @@ void Solver::run() {
             states_last_two_cycles.begin(),
             states_last_two_cycles.begin() + simparams.sim_pts_per_cycle - 1,
             states_last_two_cycles.end());
+
         last_two_cycles_time_pt_counter = simparams.sim_pts_per_cycle;
         for (size_t i = 1; i < simparams.sim_pts_per_cycle; i++) {
           state = integrator.step(state, time);
+
           states_last_two_cycles[last_two_cycles_time_pt_counter] = state;
           last_two_cycles_time_pt_counter += 1;
           interval_counter += 1;
@@ -184,8 +197,10 @@ void Solver::run() {
           }
         }
         extra_num_cycles++;
+
         converged = check_vessel_cap_convergence(states_last_two_cycles,
                                                  vessel_caps_dof_indices);
+
         assert(last_two_cycles_time_pt_counter == num_time_pts_in_two_cycles);
       }
       std::cout << "Ran simulation for " << extra_num_cycles
@@ -196,6 +211,7 @@ void Solver::run() {
         std::pair<double, double> cycle_to_cycle_errors_in_flow_and_pressure =
             get_cycle_to_cycle_errors_in_flow_and_pressure(
                 states_last_two_cycles, dof_indices);
+
         double cycle_to_cycle_error_flow =
             cycle_to_cycle_errors_in_flow_and_pressure.first;
         double cycle_to_cycle_error_pressure =
@@ -223,6 +239,13 @@ void Solver::run() {
       time -= start_time;
     }
   }
+  DEBUG_MSG("Ran time integration");
+}
+
+void Solver::run() {
+  setup_initial();
+  setup_integrator();
+  run_integration();
 }
 
 std::vector<std::pair<int, int>> Solver::get_vessel_caps_dof_indices() {
diff --git a/src/solve/Solver.h b/src/solve/Solver.h
index 4f58df727..5984a7f1d 100644
--- a/src/solve/Solver.h
+++ b/src/solve/Solver.h
@@ -32,9 +32,25 @@ class Solver {
   Solver(const nlohmann::json& config);
 
   /**
-   * @brief Run the simulation
+   * @brief Set up and initialize the simulation parameters and model
+   */
+  void setup_initial();
+
+  /**
+   * @brief Set up integrator
+   *
+   */
+  void setup_integrator();
+
+  /**
+   * @brief Run the integration
    *
    */
+  void run_integration();
+
+  /**
+   * @brief Run the simulation
+   */
   void run();
 
   /**
@@ -96,8 +112,11 @@ class Solver {
   std::shared_ptr<Model> model;
   SimulationParameters simparams;
   std::vector<State> states;
-  std::vector<double> times;
   State initial_state;
+  State state;
+  std::vector<double> times;
+  double time;
+  Integrator integrator;
 
   void sanity_checks();