diff --git a/cpp/models/ode_secir/model.h b/cpp/models/ode_secir/model.h
index dedec08aec..4f6e2749f9 100644
--- a/cpp/models/ode_secir/model.h
+++ b/cpp/models/ode_secir/model.h
@@ -24,11 +24,21 @@
 #include "memilio/compartments/simulation.h"
 #include "memilio/config.h"
 #include "memilio/epidemiology/populations.h"
+#include "memilio/io/io.h"
 #include "ode_secir/infection_state.h"
 #include "memilio/math/interpolation.h"
 #include "ode_secir/parameters.h"
 #include "memilio/math/smoother.h"
 #include "memilio/math/eigen_util.h"
+#include "unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h"
+#include <boost/math/tools/roots.hpp>
+#include <boost/math/special_functions/next.hpp> // For float_distance.
+#include <functional>
+#include <limits>
+#include <tuple> // for std::tuple and std::make_tuple.
+#include <boost/math/special_functions/cbrt.hpp> // For boost::math::cbrt.
+#include <boost/math/quadrature/gauss_kronrod.hpp>
+#include <boost/math/differentiation/finite_difference.hpp>
 
 namespace mio
 {
@@ -317,6 +327,7 @@ double get_infections_relative(const Simulation<Base>& sim, double /*t*/, const
 *@tparam Base simulation type that uses a secir compartment model. see Simulation.
 *@returns The computed reproduction number at the provided index time.
 */
+
 template <class Base>
 IOResult<ScalarType> get_reproduction_number(size_t t_idx, const Simulation<Base>& sim)
 {
@@ -522,11 +533,250 @@ IOResult<ScalarType> get_reproduction_number(size_t t_idx, const Simulation<Base
     return mio::success(eigen_vals_abs.maxCoeff());
 }
 
+IOResult<ScalarType> calc_for_lambda(Eigen::MatrixXd F, Eigen::MatrixXd V, ScalarType lambda)
+{
+    //return rho(W(t,0,lambda))
+    Eigen::MatrixXd M = 365.0 * (-V + (1 / lambda) * F);
+    //compute evolution operator
+    Eigen::MatrixXd evolution_op = M.exp();
+    //Compute the largest eigenvalue in absolute value
+    Eigen::ComplexEigenSolver<Eigen::MatrixXd> ces;
+
+    ces.compute(M);
+    const Eigen::VectorXcd eigen_vals = ces.eigenvalues();
+
+    Eigen::VectorXd eigen_vals_abs;
+    eigen_vals_abs.resize(eigen_vals.size());
+
+    for (int i = 0; i < eigen_vals.size(); i++) {
+        eigen_vals_abs[i] = std::abs(eigen_vals[i]);
+    }
+    return mio::success(eigen_vals_abs.maxCoeff());
+}
+
+template <class Base>
+IOResult<ScalarType> get_reproduction_number_with_seasonality(size_t t_idx, const Simulation<Base>& sim)
+{
+
+    if (!(t_idx < static_cast<size_t>(sim.get_result().get_num_time_points()))) {
+        return mio::failure(mio::StatusCode::OutOfRange, "t_idx is not a valid index for the TimeSeries");
+    }
+
+    auto const& params                 = sim.get_model().parameters;
+    const size_t num_groups            = (size_t)sim.get_model().parameters.get_num_groups();
+    size_t num_infected_compartments   = 5;
+    size_t total_infected_compartments = num_infected_compartments * num_groups;
+
+    Eigen::MatrixXd F(total_infected_compartments, total_infected_compartments);
+    Eigen::MatrixXd V(total_infected_compartments, total_infected_compartments);
+    F = Eigen::MatrixXd::Zero(total_infected_compartments,
+                              total_infected_compartments); //Initialize matrices F and V with zeroes
+    V = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments);
+
+    auto test_and_trace_required = 0.0;
+    auto icu_occupancy           = 0.0;
+    for (auto i = AgeGroup(0); i < (mio::AgeGroup)num_groups; ++i) {
+        auto rateINS = 0.5 / (params.template get<IncubationTime>()[i] - params.template get<SerialInterval>()[i]);
+        test_and_trace_required +=
+            (1 - params.template get<RecoveredPerInfectedNoSymptoms>()[i]) * rateINS *
+            sim.get_result().get_value(
+                t_idx)[sim.get_model().populations.get_flat_index({i, InfectionState::InfectedNoSymptoms})];
+        icu_occupancy += sim.get_result().get_value(
+            t_idx)[sim.get_model().populations.get_flat_index({i, InfectionState::InfectedCritical})];
+    }
+
+    ContactMatrixGroup const& contact_matrix = sim.get_model().parameters.template get<ContactPatterns>();
+
+    Eigen::MatrixXd cont_freq_eff(num_groups, num_groups);
+    Eigen::VectorXd divN(num_groups);
+    Eigen::VectorXd riskFromInfectedSymptomatic(num_groups);
+    Eigen::MatrixXd riskFromInfectedSymptomatic_derivatives(num_groups, num_groups);
+    Eigen::VectorXd rateINS(num_groups);
+
+    for (mio::AgeGroup k = 0; k < (mio::AgeGroup)num_groups; k++) {
+        double temp = sim.get_result().get_value(
+                          t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::Susceptible})] +
+                      sim.get_result().get_value(
+                          t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::Exposed})] +
+                      sim.get_result().get_value(
+                          t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedNoSymptoms})] +
+                      sim.get_result().get_value(
+                          t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedSymptoms})] +
+                      sim.get_result().get_value(
+                          t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedSevere})] +
+                      sim.get_result().get_value(
+                          t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedCritical})] +
+                      sim.get_result().get_value(
+                          t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::Recovered})];
+        if (temp == 0) {
+            temp = 1;
+        }
+        divN[(size_t)k] = 1 / temp;
+
+        riskFromInfectedSymptomatic[(size_t)k] = smoother_cosine(
+            test_and_trace_required, params.template get<TestAndTraceCapacity>(),
+            (params.template get<TestAndTraceCapacity>()) * 5, params.template get<RiskOfInfectionFromSymptomatic>()[k],
+            params.template get<MaxRiskOfInfectionFromSymptomatic>()[k]);
+
+        rateINS[(size_t)k] =
+            0.5 / (params.template get<IncubationTime>()[k] - params.template get<SerialInterval>()[(mio::AgeGroup)k]);
+
+        for (mio::AgeGroup l = 0; l < (mio::AgeGroup)num_groups; l++) {
+            if (test_and_trace_required < params.template get<TestAndTraceCapacity>() ||
+                test_and_trace_required > 5 * params.template get<TestAndTraceCapacity>()) {
+                riskFromInfectedSymptomatic_derivatives((size_t)k, (size_t)l) = 0;
+            }
+            else {
+                riskFromInfectedSymptomatic_derivatives((size_t)k, (size_t)l) =
+                    -0.5 * 3.14159265358979323846 *
+                    (params.template get<MaxRiskOfInfectionFromSymptomatic>()[k] -
+                     params.template get<RiskOfInfectionFromSymptomatic>()[k]) /
+                    (4 * params.template get<TestAndTraceCapacity>()) *
+                    (1 - params.template get<RecoveredPerInfectedNoSymptoms>()[l]) * rateINS[(size_t)l] *
+                    std::sin(3.14159265358979323846 / (4 * params.template get<TestAndTraceCapacity>()) *
+                             (test_and_trace_required - params.template get<TestAndTraceCapacity>()));
+            }
+        }
+
+        for (Eigen::Index l = 0; l < (Eigen::Index)num_groups; l++) {
+            cont_freq_eff(l, (size_t)k) =
+                contact_matrix.get_matrix_at(static_cast<double>(t_idx))( //Discard season_val on purpose
+                    static_cast<Eigen::Index>((size_t)l), static_cast<Eigen::Index>((size_t)k));
+        }
+    }
+
+    //Initialize the matrix F
+    for (size_t i = 0; i < num_groups; i++) {
+
+        for (size_t j = 0; j < num_groups; j++) {
+
+            double temp = 0;
+            for (Eigen::Index k = 0; k < (Eigen::Index)num_groups; k++) {
+                temp += cont_freq_eff(i, k) *
+                        sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
+                            {(mio::AgeGroup)k, InfectionState::InfectedSymptoms})] *
+                        riskFromInfectedSymptomatic_derivatives(k, j) * divN[k];
+            }
+
+            F(i, j + num_groups) =
+                sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
+                    {(mio::AgeGroup)i, InfectionState::Susceptible})] *
+                params.template get<TransmissionProbabilityOnContact>()[(mio::AgeGroup)i] *
+                (cont_freq_eff(i, j) * params.template get<RelativeTransmissionNoSymptoms>()[(mio::AgeGroup)j] *
+                     divN[(size_t)j] +
+                 temp);
+        }
+
+        for (size_t j = 0; j < num_groups; j++) {
+            F(i, j + 2 * num_groups) = sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
+                                           {(mio::AgeGroup)i, InfectionState::Susceptible})] *
+                                       params.template get<TransmissionProbabilityOnContact>()[(mio::AgeGroup)i] *
+                                       cont_freq_eff(i, j) * riskFromInfectedSymptomatic[(size_t)j] * divN[(size_t)j];
+        }
+    }
+
+    //Initialize the matrix V
+    for (Eigen::Index i = 0; i < (Eigen::Index)num_groups; i++) {
+
+        double rateE = 1.0 / (2 * params.template get<SerialInterval>()[(mio::AgeGroup)i] -
+                              params.template get<IncubationTime>()[(mio::AgeGroup)i]);
+
+        double criticalPerSevereAdjusted = smoother_cosine(
+            icu_occupancy, 0.90 * params.template get<ICUCapacity>(), params.template get<ICUCapacity>(),
+            params.template get<CriticalPerSevere>()[(mio::AgeGroup)i], 0);
+
+        V(i, i)                           = rateE;
+        V(i + num_groups, i)              = -rateE;
+        V(i + num_groups, i + num_groups) = rateINS[i];
+        V(i + 2 * num_groups, i + num_groups) =
+            -(1 - params.template get<RecoveredPerInfectedNoSymptoms>()[(mio::AgeGroup)i]) * rateINS[i];
+        V(i + 2 * num_groups, i + 2 * num_groups) = (1 / params.template get<TimeInfectedSymptoms>()[(mio::AgeGroup)i]);
+        V(i + 3 * num_groups, i + 2 * num_groups) =
+            -params.template get<SeverePerInfectedSymptoms>()[(mio::AgeGroup)i] /
+            params.template get<TimeInfectedSymptoms>()[(mio::AgeGroup)i];
+        V(i + 3 * num_groups, i + 3 * num_groups) = 1 / (params.template get<TimeInfectedSevere>()[(mio::AgeGroup)i]);
+        V(i + 4 * num_groups, i + 3 * num_groups) =
+            -criticalPerSevereAdjusted / (params.template get<TimeInfectedSevere>()[(mio::AgeGroup)i]);
+
+        V(i + 4 * num_groups, i + 4 * num_groups) = 1 / (params.template get<TimeInfectedCritical>()[(mio::AgeGroup)i]);
+
+        if (!(icu_occupancy < 0.9 * params.template get<ICUCapacity>() ||
+              icu_occupancy > (double)(params.template get<ICUCapacity>()))) {
+            for (size_t j = 0; j < num_groups; j++) {
+                V(i + 4 * num_groups, j + 4 * num_groups) -=
+                    sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
+                        {(mio::AgeGroup)i, InfectionState::InfectedSevere})] /
+                    params.template get<TimeInfectedSevere>()[(mio::AgeGroup)i] * 5 *
+                    params.template get<CriticalPerSevere>()[(mio::AgeGroup)i] * 3.141592653589793 /
+                    (params.template get<ICUCapacity>()) *
+                    std::sin(3.141592653589793 / (0.1 * params.template get<ICUCapacity>()) *
+                             (icu_occupancy - 0.9 * params.template get<ICUCapacity>()));
+            }
+        }
+    }
+
+    //Set up the iteration
+    ScalarType start_guess                                  = get_reproduction_number(t_idx, sim);
+    boost::uintmax_t maxit                                  = 30;
+    std::function<std::vector<double>(ScalarType)> calc_val = [F, V](ScalarType lambda) {
+        //return spectral radius and derivative of spectral radius
+        //return rho(W(t,0,lambda))
+        Eigen::MatrixXd M = 365.0 * (-V + (1 / lambda) * F);
+        //compute evolution operator
+        Eigen::MatrixXd evolution_op = M.exp();
+        //Compute the largest eigenvalue in absolute value
+        Eigen::ComplexEigenSolver<Eigen::MatrixXd> ces;
+
+        ces.compute(M);
+        const Eigen::VectorXcd eigen_vals = ces.eigenvalues();
+        const Eigen::MatrixXcd eigen_vecs = ces.eigenvectors();
+
+        Eigen::VectorXcd eigen_vals_abs;
+        eigen_vals_abs.resize(eigen_vals.size());
+
+        for (int i = 0; i < eigen_vals.size(); i++) {
+            eigen_vals_abs[i] = std::abs(eigen_vals[i]);
+        }
+        auto it_max_ev =
+            std::find(eigen_vals_abs.data(), eigen_vals_abs.data() + eigen_vals_abs.size(), eigen_vals_abs.maxCoeff());
+
+        auto idx_max_ev                = std::distance(eigen_vals_abs.data(), it_max_ev);
+        Eigen::VectorXcd max_eigen_vec = eigen_vecs.col(idx_max_ev);
+
+        //Now calculate derivative of exp(A(t)) evaluated at max_eigen_vec at lambda
+        Eigen::MatrixXcd MatDeriv(F.rows(), F.cols());
+        for (Eigen::Index i = 0; i < F.rows(); i++) {
+            for (Eigen::Index j = 0; j < F.cols(); j++) {
+                std::function<std::complex<double>(std::complex<double>)> integrand = [F, evolution_op, lambda, i,
+                                                                                       j](ScalarType alpha) {
+                    Eigen::MatrixXcd result(F.rows(), F.cols());
+                    result = -(1 / (lambda * lambda)) * (alpha * evolution_op).exp() * F *
+                             ((1 - alpha) * evolution_op).exp();
+                    return result(i, j);
+                };
+                MatDeriv(i, j) =
+                    boost::math::quadrature::gauss_kronrod<std::complex<ScalarType>, 100>::integrate(integrand, 0, 1);
+            }
+        }
+
+        Eigen::VectorXcd term1(F.rows());
+        term1 = MatDeriv * max_eigen_vec;
+
+        //Calculate derivative of largest eigenvector numerically
+
+        Eigen::VectorXcd dx_max_eigen_vec(F.rows());
+        Eigen::VectorXd term2 = evolution_op * dx_max_eigen_vec;
+
+        std::function<ScalarType(ScalarType)> deriv_max_eigen_vec = return std::make_tuple(*(it_max_ev), );
+    };
+}
+
 /**
 *@brief Computes the reproduction number for all time points of the Model output obtained by the Simulation.
 *@param sim The Model Simulation.
 *@returns vector containing all reproduction numbers
 */
+
 template <class Base>
 Eigen::VectorXd get_reproduction_numbers(const Simulation<Base>& sim)
 {