diff --git a/ComponentLib/CMakeLists.txt b/ComponentLib/CMakeLists.txt
index a97c973..16b5a24 100644
--- a/ComponentLib/CMakeLists.txt
+++ b/ComponentLib/CMakeLists.txt
@@ -62,4 +62,5 @@
 
 add_subdirectory(PowerFlow)
 add_subdirectory(PowerElectronics)
+add_subdirectory(DynamicPhasor)
 
diff --git a/ComponentLib/DynamicPhasor/Bus/CMakeLists.txt b/ComponentLib/DynamicPhasor/Bus/CMakeLists.txt
new file mode 100644
index 0000000..8d713c9
--- /dev/null
+++ b/ComponentLib/DynamicPhasor/Bus/CMakeLists.txt
@@ -0,0 +1,5 @@
+gridkit_add_library(slackbus
+  SOURCES
+    SlackBus.cpp
+  OUTPUT_NAME
+    gridkit_slackbus)
\ No newline at end of file
diff --git a/ComponentLib/DynamicPhasor/Bus/PhasorBus.hpp b/ComponentLib/DynamicPhasor/Bus/PhasorBus.hpp
new file mode 100644
index 0000000..3e8ac00
--- /dev/null
+++ b/ComponentLib/DynamicPhasor/Bus/PhasorBus.hpp
@@ -0,0 +1,100 @@
+/* Paul Moon 6/26/2024 */
+
+#ifndef _PHASOR_BUS_HPP_
+#define _PHASOR_BUS_HPP_
+
+#include <ModelEvaluatorImpl.hpp>
+
+namespace ModelLib
+{
+    /*!
+     * @brief Base class for all dynamic phasor buses.
+     *
+     * Derived bus types:
+     *   0 - swing bus (Vr and Vi are constants)
+     *   1 - PV bus    (P and V are constants)
+     *   2 - PQ bus    (P and Q are constants)
+     *
+     * @todo Consider static instead of dynamic polymorphism for
+     * bus types. Create Bus class that takes template parameter
+     * BusType.
+     */
+    template  <class ScalarT, typename IdxT>
+    class PhasorBus : public ModelEvaluatorImpl<ScalarT, IdxT>
+    {
+    protected:
+        using ModelEvaluatorImpl<ScalarT, IdxT>::size_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::nnz_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::time_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::alpha_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::rtol_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::atol_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::y_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::yp_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::tag_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::f_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::g_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::yB_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::ypB_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::fB_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::gB_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::param_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::param_up_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::param_lo_;
+
+    public:
+        typedef typename ModelEvaluatorImpl<ScalarT, IdxT>::real_type real_type;
+
+        enum BusType{PQ=1, PV, Slack, Isolated};
+
+        PhasorBus(IdxT id) : busID_(id) {}
+        virtual ~PhasorBus(){}
+
+        // Set defaults for ModelEvaluator methods
+        virtual int allocate() { return 0;}
+        virtual int initialize() { return 0;}
+        virtual int tagDifferentiable() { return 0;}
+        virtual int evaluateResidual() { return 0;}
+        virtual int evaluateJacobian() { return 0;}
+        virtual int evaluateIntegrand() { return 0;}
+
+        virtual int initializeAdjoint() { return 0;}
+        virtual int evaluateAdjointResidual() { return 0;}
+        //virtual int evaluateAdjointJacobian() { return 0;}
+        virtual int evaluateAdjointIntegrand() { return 0;}
+        virtual void updateTime(real_type, real_type) {} // <- throw exception here
+
+        // Pure virtual methods specific to Bus types
+        virtual ScalarT& Vr() = 0;
+        virtual const ScalarT& Vr() const = 0;
+        virtual ScalarT& Vi() = 0;
+        virtual const ScalarT& Vi() const = 0;
+        virtual ScalarT& P() = 0;
+        virtual const ScalarT& P() const = 0;
+        virtual ScalarT& Q() = 0;
+        virtual const ScalarT& Q() const = 0;
+
+        virtual ScalarT& lambdaP() = 0;
+        virtual const ScalarT& lambdaP() const = 0;
+        virtual ScalarT& lambdaQ() = 0;
+        virtual const ScalarT& lambdaQ() const = 0;
+        virtual ScalarT& PB() = 0;
+        virtual const ScalarT& PB() const = 0;
+        virtual ScalarT& QB() = 0;
+        virtual const ScalarT& QB() const = 0;
+
+        virtual const int BusType() const = 0;
+
+        virtual const IdxT BusID() const
+        {
+            return busID_;
+        }
+
+    protected:
+        const IdxT busID_;
+    }; // class PhasorBus
+
+} // namespace ModelLib
+
+
+#endif // _PHASOR_BUS_HPP_
diff --git a/ComponentLib/DynamicPhasor/Bus/SlackBus.cpp b/ComponentLib/DynamicPhasor/Bus/SlackBus.cpp
new file mode 100644
index 0000000..d4924aa
--- /dev/null
+++ b/ComponentLib/DynamicPhasor/Bus/SlackBus.cpp
@@ -0,0 +1,87 @@
+/* Created by Paul Moon 6/26/2024 */
+#include <iostream>
+#include <cmath>
+#include "SlackBus.hpp"
+
+namespace ModelLib {
+
+/*!
+ * @brief Constructor for a slack bus
+ *
+ * Arguments passed to ModelEvaluatorImpl:
+ * - Number of equations = 0 (size_)
+ * - Number of variables = 0 (size_)
+ * - Number of quadratures = 0
+ * - Number of optimization parameters = 0
+ */
+template <class ScalarT, typename IdxT>
+SlackBus<ScalarT, IdxT>::SlackBus()
+    : PhasorBus<ScalarT, IdxT>(0), Vr_(0.0), Vi_(0.0), P_(0.0), Q_(0.0), PB_(0.0), QB_(0.0)
+{
+    //std::cout << "Create SlackBus..." << std::endl;
+    //std::cout << "Number of equations is " << size_ << std::endl;
+
+    size_ = 0;
+}
+
+/*!
+ * @brief SlackBus constructor.
+ *
+ * Arguments passed to ModelEvaluatorImpl:
+ * - Number of equations = 0 (size_)
+ * - Number of variables = 0 (size_)
+ * - Number of quadratures = 0
+ * - Number of optimization parameters = 0
+ */
+template <class ScalarT, typename IdxT>
+SlackBus<ScalarT, IdxT>::SlackBus(ScalarT Vr, ScalarT Vi)
+    : PhasorBus<ScalarT, IdxT>(0), Vr_(Vr), Vi_(Vi), P_(0.0), Q_(0.0), PB_(0.0), QB_(0.0)
+{
+    //std::cout << "Create SlackBus..." << std::endl;
+    //std::cout << "Number of equations is " << size_ << std::endl;
+    P() = 0.0;
+    Q() = 0.0;
+    size_ = 0;
+}
+
+template <class ScalarT, typename IdxT>
+SlackBus<ScalarT, IdxT>::SlackBus(BusData& data)
+    : PhasorBus<ScalarT, IdxT>(data.bus_i), Vr_(data.Vm*cos(data.Va)), Vi_(data.Vm*sin(data.Va))
+{
+    //std::cout << "Create SlackBus..." << std::endl;
+    //std::cout << "Number of equations is " << size_ << std::endl;
+    P() = 0.0;
+    Q() = 0.0;
+    size_ = 0;
+}
+
+template <class ScalarT, typename IdxT>
+SlackBus<ScalarT, IdxT>::~SlackBus()
+{
+}
+
+template <class ScalarT, typename IdxT>
+int SlackBus<ScalarT, IdxT>::evaluateResidual()
+{
+    // std::cout << "Evaluating residual of a slack bus ...\n";
+    P() = 0.0;
+    Q() = 0.0;
+    return 0;
+}
+
+template <class ScalarT, typename IdxT>
+int SlackBus<ScalarT, IdxT>::evaluateAdjointResidual()
+{
+    PB() = 0.0;
+    QB() = 0.0;
+    return 0;
+}
+
+
+// Available template instantiations
+template class SlackBus<double, long int>;
+template class SlackBus<double, size_t>;
+
+
+} // namespace ModelLib
+
diff --git a/ComponentLib/DynamicPhasor/Bus/SlackBus.hpp b/ComponentLib/DynamicPhasor/Bus/SlackBus.hpp
new file mode 100644
index 0000000..c6bb630
--- /dev/null
+++ b/ComponentLib/DynamicPhasor/Bus/SlackBus.hpp
@@ -0,0 +1,146 @@
+/* Created by Paul Moon 6/26/2024 */
+#ifndef _SLACK_BUS_HPP_
+#define _SLACK_BUS_HPP_
+
+#include <PowerSystemData.hpp>
+#include "PhasorBus.hpp"
+
+namespace ModelLib
+{
+    /*!
+     * @brief Implementation of a slack bus.
+     *
+     * Slack bus sets voltage _V_ and phase _theta_ as constants.
+     * Active and reactive power, _P_ and _Q_, are component model outputs,
+     * but are computed outside the SlackBus class.
+     *
+     *
+     */
+    template  <class ScalarT, typename IdxT>
+    class SlackBus : public PhasorBus<ScalarT, IdxT>
+    {
+        using PhasorBus<ScalarT, IdxT>::size_;
+        using PhasorBus<ScalarT, IdxT>::y_;
+        using PhasorBus<ScalarT, IdxT>::yp_;
+        using PhasorBus<ScalarT, IdxT>::f_;
+        using PhasorBus<ScalarT, IdxT>::g_;
+        using PhasorBus<ScalarT, IdxT>::atol_;
+        using PhasorBus<ScalarT, IdxT>::rtol_;
+
+    public:
+        using real_type = typename ModelEvaluatorImpl<ScalarT, IdxT>::real_type;
+        using BusData = GridKit::PowerSystemData::BusData<real_type, IdxT>;
+
+        SlackBus();
+        SlackBus(ScalarT Vr, ScalarT Vi);
+        SlackBus(BusData& data);
+        virtual ~SlackBus();
+        virtual int evaluateResidual();
+        virtual int evaluateAdjointResidual();
+
+        /// @todo Should slack bus allow changing voltage?
+        virtual ScalarT& Vr()
+        {
+            return Vr_;
+        }
+
+        virtual const ScalarT& Vr() const
+        {
+            return Vr_;
+        }
+
+        /// @todo Should slack bus allow changing phase?
+        virtual ScalarT& Vi()
+        {
+            return Vi_;
+        }
+
+        virtual const ScalarT& Vi() const
+        {
+            return Vi_;
+        }
+
+        virtual ScalarT& P()
+        {
+            return P_;
+        }
+
+        virtual const ScalarT& P() const
+        {
+            return P_;
+        }
+
+        virtual ScalarT& Q()
+        {
+            return Q_;
+        }
+
+        virtual const ScalarT& Q() const
+        {
+            return Q_;
+        }
+
+        /// @todo Should slack bus allow changing voltage?
+        virtual ScalarT& lambdaP()
+        {
+            return thetaB_;
+        }
+
+        virtual const ScalarT& lambdaP() const
+        {
+            return thetaB_;
+        }
+
+        /// @todo Should slack bus allow changing phase?
+        virtual ScalarT& lambdaQ()
+        {
+            return VB_;
+        }
+
+        virtual const ScalarT& lambdaQ() const
+        {
+            return VB_;
+        }
+
+        virtual ScalarT& PB()
+        {
+            return PB_;
+        }
+
+        virtual const ScalarT& PB() const
+        {
+            return PB_;
+        }
+
+        virtual ScalarT& QB()
+        {
+            return QB_;
+        }
+
+        virtual const ScalarT& QB() const
+        {
+            return QB_;
+        }
+
+        virtual const int BusType() const
+        {
+            return PhasorBus<ScalarT, IdxT>::BusType::Slack;
+        }
+
+    private:
+        ScalarT Vr_;
+        ScalarT Vi_;
+        ScalarT P_;
+        ScalarT Q_;
+
+        ScalarT VB_;
+        ScalarT thetaB_;
+        ScalarT PB_;
+        ScalarT QB_;
+
+    }; // class SlackBus
+
+} // namespace ModelLib
+
+
+#endif // _SLACK_BUS_HPP_
\ No newline at end of file
diff --git a/ComponentLib/DynamicPhasor/CMakeLists.txt b/ComponentLib/DynamicPhasor/CMakeLists.txt
new file mode 100644
index 0000000..00dc18e
--- /dev/null
+++ b/ComponentLib/DynamicPhasor/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_subdirectory(SynchronousMachine)
+add_subdirectory(Governor)
+add_subdirectory(Bus)
\ No newline at end of file
diff --git a/ComponentLib/DynamicPhasor/Governor/CMakeLists.txt b/ComponentLib/DynamicPhasor/Governor/CMakeLists.txt
new file mode 100644
index 0000000..047c0d2
--- /dev/null
+++ b/ComponentLib/DynamicPhasor/Governor/CMakeLists.txt
@@ -0,0 +1,5 @@
+gridkit_add_library(governor
+  SOURCES
+    Governor.cpp
+  OUTPUT_NAME
+    gridkit_governor)
\ No newline at end of file
diff --git a/ComponentLib/DynamicPhasor/Governor/Governor.cpp b/ComponentLib/DynamicPhasor/Governor/Governor.cpp
new file mode 100644
index 0000000..42e13d8
--- /dev/null
+++ b/ComponentLib/DynamicPhasor/Governor/Governor.cpp
@@ -0,0 +1,360 @@
+// Created by Paul Moon 7/11/2024
+
+#include <iostream>
+#include <cmath>
+#include <complex>
+
+#include <ComponentLib/PowerFlow/Bus/BaseBus.hpp>
+#include "Governor.hpp"
+
+namespace ModelLib
+{
+
+    /*!
+     * @brief Constructor for a simple governor model
+     *
+     * Arguments passed to ModelEvaluatorImpl:
+     * - Number of equations = 6 differential + 9 algebraic = 15
+     * - Quadratures and optimization parameters are not being used.
+     */
+    template <class ScalarT, typename IdxT>
+    Governor<ScalarT, IdxT>::Governor(bus_type *bus, ScalarT P0, ScalarT Q0)
+        : ModelEvaluatorImpl<ScalarT, IdxT>(18, 0, 0),
+          H_(3.0),                    // Intertia constant
+          D_(0.0),                    // Damping coefficient
+          Xq_(0.5),                   // Machine reactance parameter
+          Xd_(2.1),                   // Machine reactance parameter
+          Xqp_(0.5),                  // Machine reactance parameter
+          Xdp_(0.2),                  // Machine reactance parameter
+          Xqpp_(0.18),                // Machine reactance parameter
+          Xdpp_(0.18),                // Machine reactance parameter
+          Rs_(0.0),                   // Winding resistance
+          Tq0p_(0.75),                // Time constant
+          Td0p_(7.0),                 // Time constant
+          Tq0pp_(0.05),               // Time constant
+          Td0pp_(0.035),              // Time constant
+          Ef_(1.8336),                // Field winding voltage from the excitation system. An exciter would set this value
+          X_(0.22),                   // Machine reactance parameter
+          Xl_(0.15),                  // Machine reactance parameter
+          omega_s_(0.0),              // Relative speed deviation. 0 is synchronous, -1 is 0 Hz.
+          omega0_(2.0 * 60.0 * M_PI), // Nominal frequency
+          S1_(0.05),                  // Saturation parameter
+          S12_(0.2),                  // Saturation parameter
+          R_(0.05),  // Droop constant
+          T1_(0.5),  // Valve time delay
+          T2_(2.5),  // Turbine numerator time constant
+          T3_(7.5),  // Turbine delay
+          Vmax_(1),  // Max valve position
+          Vmin_(-1), // Min valve position
+          Dt_(0),    // Turbine damping
+          P0_(P0),  // Real power
+          Q0_(Q0),  // Reactive power
+          offsetGov_(15),
+          bus_(bus)
+    {
+    }
+
+    template <class ScalarT, typename IdxT>
+    Governor<ScalarT, IdxT>::~Governor()
+    {
+    }
+
+    /*!
+     * @brief This function will be used to allocate sparse Jacobian matrices.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int Governor<ScalarT, IdxT>::allocate()
+    {
+        tag_.resize(size_);
+
+        return 0;
+    }
+
+    /**
+     * @brief Initialization of the governor model
+     *
+     * Initialization equations are derived from Adam
+     * Birchfield's paper on power electronics.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int Governor<ScalarT, IdxT>::initialize()
+    {
+        // Compute values for saturation constants and real/imaginary voltage values
+        double SB = ((1 / 0.2) * (1 / 0.2)) * ((sqrt(S1_) - sqrt(S12_)) * (sqrt(S1_) - sqrt(S12_)));
+        double SA = 1 - sqrt(S1_ / SB);
+        const double Vr = V() * cos(theta());
+        const double Vi = V() * sin(theta());
+
+        const double Ir = P0_ / V();
+        const double Ii = 0.0;
+
+        // Compute initial guess for the generator voltage phase
+        const double delta = atan((Xq_ * Ir + Vi) / (-Xq_ * Ii + Vr));
+
+        // Compute initial guesses for generator currents and potentials in d-q frame
+        const double Id = Ir * sin(delta) - Ii * cos(delta);
+        const double Iq = Ir * cos(delta) + Ii * sin(delta);
+
+        const double Vq = Vr * cos(delta) + Vi * sin(delta) + Id * Xqpp_ + Iq * Rs_;
+        const double Vd = Vr * sin(delta) - Vi * cos(delta) + Id * Rs_ - Iq * Xqpp_;
+
+        // Compute initial guesses for flux
+        const double Psiqpp = -Vd / (1 + omega_s_);
+        const double Psidpp = Vq / (1 + omega_s_);
+
+        const double Psipp = sqrt((Psiqpp * Psiqpp) + (Psidpp * Psidpp));
+        const double ksat = SB * ((Psipp - SA) * (Psipp - SA));
+
+        // Compute initial guess for Telec. Note that Pmech would normally be changed by a governor.
+        const double Telec = (Psidpp - Id * Xdpp_) * Iq - (Psiqpp - Iq * Xdpp_) * Id;
+        const double Pmech = Telec;
+
+        // Compute remaining initial guesses.
+        const double Psiqp = (1.0 / (1.0 + (Xqp_ - Xqpp_) / (Xqpp_ - Xl_))) * (-Psiqpp * (Xqp_ - Xl_) / (Xqpp_ - Xl_) + Iq * (Xqp_ - Xl_));
+
+        const double Edp = Psiqp - Iq * (Xqp_ - Xl_);
+
+        const double Psidp = -(1.0 / (1.0 + (Xdp_ - Xdpp_) / (Xdpp_ - Xl_))) * (Id * (Xdp_ - Xl_) - Psidpp * (Xdp_ - Xl_) / (Xdpp_ - Xl_));
+
+        const double Eqp = Psidp + Id * (Xdp_ - Xl_);
+
+        const double Pv_init = Vmin_;
+        const double Ptx_init = Pmech * T3_ - T2_ * Pv_init-T3_*Dt_*omega_s_;
+        const double Pref_init = P0_;
+
+        y_[0] = delta;
+        y_[1] = omega_s_;
+        y_[2] = Eqp;
+        y_[3] = Psidp;
+        y_[4] = Psiqp;
+        y_[5] = Edp;
+        y_[6] = Psiqpp;
+        y_[7] = Psidpp;
+        y_[8] = Psipp;
+        y_[9] = ksat;
+        y_[10] = Vd;
+        y_[11] = Vq;
+        y_[12] = Telec;
+        y_[13] = Id;
+        y_[14] = Iq;
+        y_[15] = Pmech;
+
+        yp_[0] = 0.0;
+        yp_[1] = 0.0;
+        yp_[2] = 0.0;
+        yp_[3] = 0.0;
+        yp_[4] = 0.0;
+        yp_[5] = 0.0;
+        yp_[6] = 0.0;
+        yp_[7] = 0.0;
+        yp_[8] = 0.0;
+        yp_[9] = 0.0;
+        yp_[10] = 0.0;
+        yp_[11] = 0.0;
+        yp_[12] = 0.0;
+        yp_[13] = 0.0;
+        yp_[14] = 0.0;
+        yp_[15] = 0.0;
+
+        y_[offsetGov_ + 1] = Ptx_init;
+        y_[offsetGov_ + 2] = Pv_init;
+
+        yp_[offsetGov_ + 1] = 0.0;
+        yp_[offsetGov_ + 2] = 0.0;
+
+        for (int i = 0; i < size_; ++i) {
+            std::cout << y_[i] << ", ";
+        }
+        std::cout << std::endl;
+
+        return 0;
+    }
+
+    /**
+     * \brief Identify differential variables.
+     */
+    template <class ScalarT, typename IdxT>
+    int Governor<ScalarT, IdxT>::tagDifferentiable()
+    {
+        tag_[0] = true;
+        tag_[1] = true;
+        tag_[2] = true;
+        tag_[3] = true;
+        tag_[4] = true;
+        tag_[5] = true;
+        for (IdxT i = 6; i < size_; ++i)
+        {
+            tag_[i] = false;
+        }
+        tag_[offsetGov_ + 1] = true;
+        tag_[offsetGov_ + 2] = true;
+
+        return 0;
+    }
+
+    /**
+     * @brief Computes residual vector for the generator model.
+     *
+     * Residual equations are given by Adam Birchfield's paper
+     * on power electronics.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int Governor<ScalarT, IdxT>::evaluateResidual()
+    {
+        // These constants are the same as in the initialization state. Maybe move them to .hpp file?
+        double SB = ((1 / 0.2) * (1 / 0.2)) * ((sqrt(S1_) - sqrt(S12_)) * (sqrt(S1_) - sqrt(S12_)));
+        double SA = 1 - sqrt(S1_ / SB);
+        const double Vr = V() * cos(theta());
+        const double Vi = V() * sin(theta());
+        double Pref = 1.0;
+
+        // The first 6 functions are differential equations (involving yp_)
+        f_[0] = -yp_[0] + omega0_ * y_[1];
+        f_[1] = -yp_[1] + (1 / (2 * H_)) * ((y_[15] - D_ * y_[1]) / (1 + y_[1]) - y_[12]);
+        f_[2] = -yp_[2] * Td0p_ + Ef_ - (y_[2] + (Xd_ - Xdp_) * (y_[13] + ((Xdp_ - Xdpp_) / ((Xdp_ - Xl_) * (Xdp_ - Xl_))) * (y_[2] - y_[3] - (Xdp_ - Xl_) * y_[13])) + y_[7] * y_[9]);
+        f_[3] = -yp_[3] + (1 / Td0pp_) * (y_[2] - y_[3] - (Xdp_ - Xl_) * y_[13]);
+        f_[4] = -yp_[4] + (1 / Tq0pp_) * (y_[5] - y_[4] + (Xqp_ - Xl_) * y_[14]);
+        f_[5] = -yp_[5] + (1 / Tq0p_) * (-y_[5] + (((Xq_ - Xl_) / (Xd_ - Xl_)) * y_[6] * y_[9] + (Xq_ - Xqp_) * (y_[14] - ((Xqp_ - Xqpp_) / ((Xqp_ - Xl_) * (Xqp_ - Xl_)) * y_[5] + y_[14]) * (Xqp_ - Xl_) - y_[4])));
+
+        // Remaining equations are all algebraic
+        f_[6] = -y_[6] - y_[4] * ((Xqp_ - Xqpp_) / (Xqp_ - Xl_)) - y_[5] * ((Xqpp_ - Xl_) / (Xqp_ - Xl_));
+        f_[7] = -y_[7] + y_[3] * ((Xdp_ - Xdpp_) / (Xdp_ - Xl_)) + y_[2] * ((Xdpp_ - Xl_) / (Xdp_ - Xl_));
+        f_[8] = -y_[8] + sqrt((y_[7] * y_[7]) + (y_[6] * y_[6]));
+        f_[9] = -y_[9] + SB * ((y_[8] - SA) * (y_[8] - SA));
+        f_[10] = -y_[10] - y_[6] * (1 + y_[1]);
+        f_[11] = -y_[11] + y_[7] * (1 + y_[1]);
+        f_[12] = -y_[12] + y_[14] * (y_[7] - y_[13] * Xdpp_) - y_[13] * (y_[6] - y_[14] * Xdpp_);
+        f_[13] = Vr * sin(y_[0]) - Vi * cos(y_[0]) + y_[13] * Rs_ - y_[14] * Xqpp_;
+        f_[14] = Vr * cos(y_[0]) - Vi * sin(y_[0]) + y_[13] * Xqpp_ - y_[14] * Rs_;
+        f_[15] = -y_[offsetGov_] + (1.0 / T3_) * (y_[offsetGov_ + 1] + T2_ * y_[offsetGov_ + 2]) - Dt_ * y_[1];
+
+        P() += Pg();
+        Q() += Qg();
+
+        if (y_[offsetGov_ + 2] >= Vmax_) 
+        {
+            y_[offsetGov_ + 2] = Vmax_;
+            if (((1.0/T1_) * (-y_[offsetGov_ + 2]) + (1.0/R_) * (Pref - y_[1])) > 0)
+            {
+                f_[offsetGov_+1] = -yp_[offsetGov_ + 2];
+            }
+            else
+            {
+                f_[offsetGov_+1] = -yp_[offsetGov_ + 2] + (1.0/T1_) * (-y_[offsetGov_ + 2] + (1.0/R_) * (Pref - y_[1]));
+            }
+        }
+        else if (y_[offsetGov_ + 2] <= Vmin_) 
+        {
+            y_[offsetGov_ + 2] = Vmin_;
+            if (((1.0/T1_) * (-y_[offsetGov_ + 2]) + (1.0/R_) * (Pref - y_[1])) > 0)
+            {
+                f_[offsetGov_+1] = -yp_[offsetGov_ + 2];
+            }
+            else
+            {
+                f_[offsetGov_+1] = -yp_[offsetGov_ + 2] + (1.0/T1_) * (-y_[offsetGov_ + 2] + (1.0/R_) * (Pref - y_[1]));
+            }
+        }
+        else
+        {
+            f_[offsetGov_+1] = -yp_[offsetGov_ + 2] + (1.0/T1_) * (-y_[offsetGov_ + 2] + (1.0/R_) * (Pref - y_[1]));
+        }
+        f_[offsetGov_+2] = -yp_[offsetGov_ + 1] + y_[offsetGov_ + 2] - (1.0 / T3_) * (y_[offsetGov_ + 1] + T2_ * y_[offsetGov_ + 2]);
+        //f_[offsetGov_+3] = -y_[offsetGov_] + (1.0 / T3_) * (y_[offsetGov_ + 1] + T2_ * y_[offsetGov_ + 2]) - Dt_ * y_[1];
+        return 0;
+    }
+
+    /**
+     * @brief Creates the Jacobian using J and Jp.
+     *
+     * Values were calculated by hand and verified with MATLAB.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int Governor<ScalarT, IdxT>::evaluateJacobian()
+    {
+        return 0;
+    }
+
+    template <class ScalarT, typename IdxT>
+    ScalarT Governor<ScalarT, IdxT>::Pg()
+    {
+        return y_[14]*y_[2] + y_[13]*y_[5] + (Xqp_ - Xdp_)*y_[13]*y_[14] - Rs_*(y_[13]*y_[13] + y_[14]*y_[14]);
+    }
+
+    /*y_[0] = delta;
+        y_[1] = omega_s_;
+        y_[2] = Eqp;
+        y_[3] = Psidp;
+        y_[4] = Psiqp;
+        y_[5] = Edp;
+        y_[6] = Psiqpp;
+        y_[7] = Psidpp;
+        y_[8] = Psipp;
+        y_[9] = ksat;
+        y_[10] = Vd;
+        y_[11] = Vq;
+        y_[12] = Telec;
+        y_[13] = Id;
+        y_[14] = Iq;
+        y_[15] = Pmech;*/
+
+/**
+ * Generator reactive power Qg.
+ *
+ * \f[ Q_g = E_q' I_d - E_d' I_q - X_d' I_d^2 - X_q' I_q^2 \f]
+ */
+    template <class ScalarT, typename IdxT>
+    ScalarT Governor<ScalarT, IdxT>::Qg()
+    {
+        return -y_[14]*y_[5] + y_[13]*y_[2] - Xdp_*y_[13]*y_[13] - Xqp_*y_[14]*y_[14];
+    }
+
+    /**
+     * @brief This is not used, but is implemented to keep SystemModel.hpp happy.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int Governor<ScalarT, IdxT>::evaluateIntegrand()
+    {
+        return 0;
+    }
+
+    /**
+     * @brief This is not used, but is implemented to keep SystemModel.hpp happy.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int Governor<ScalarT, IdxT>::initializeAdjoint()
+    {
+        return 0;
+    }
+
+    /**
+     * @brief This is not used, but is implemented to keep SystemModel.hpp happy.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int Governor<ScalarT, IdxT>::evaluateAdjointResidual()
+    {
+        return 0;
+    }
+
+    /**
+     * @brief This is not used, but is implemented to keep SystemModel.hpp happy.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int Governor<ScalarT, IdxT>::evaluateAdjointIntegrand()
+    {
+        return 0;
+    }
+
+    // Available template instantiations
+    template class Governor<double, long int>;
+    template class Governor<double, size_t>;
+
+} // namespace ModelLib
\ No newline at end of file
diff --git a/ComponentLib/DynamicPhasor/Governor/Governor.hpp b/ComponentLib/DynamicPhasor/Governor/Governor.hpp
new file mode 100644
index 0000000..98c0919
--- /dev/null
+++ b/ComponentLib/DynamicPhasor/Governor/Governor.hpp
@@ -0,0 +1,149 @@
+// Created by Paul Moon 7/11/2024
+
+#ifndef _GOVERNOR_H_
+#define _GOVERNOR_H_
+
+#include <ModelEvaluatorImpl.hpp>
+
+namespace ModelLib
+{
+    template <class ScalarT, typename IdxT> class BaseBus;
+}
+
+namespace ModelLib
+{
+    /*!
+     * @brief Implementation of a TGOV1 Turbine-Governor model.
+     *
+     */
+    template  <class ScalarT, typename IdxT>
+    class Governor : public ModelEvaluatorImpl<ScalarT, IdxT>
+    {
+        using ModelEvaluatorImpl<ScalarT, IdxT>::size_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::nnz_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::time_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::alpha_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::y_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::yp_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::tag_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::f_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::g_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::yB_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::ypB_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::fB_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::gB_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::param_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::param_up_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::param_lo_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::J_;
+
+        typedef typename ModelEvaluatorImpl<ScalarT, IdxT>::real_type real_type;
+        typedef BaseBus<ScalarT, IdxT> bus_type;
+
+    public:
+        Governor(BaseBus<ScalarT, IdxT>* bus, ScalarT P0, ScalarT Q0);
+        virtual ~Governor();
+
+        int allocate();
+        int initialize();
+        int tagDifferentiable();
+        int evaluateResidual();
+        int evaluateJacobian();
+
+        int evaluateIntegrand();
+        int initializeAdjoint();
+        int evaluateAdjointResidual();
+        int evaluateAdjointIntegrand();
+
+        ScalarT Pg();
+        ScalarT Qg();
+
+
+        void updateTime(real_type t, real_type a)
+        {
+            time_ = t;
+            alpha_ = a;
+        }
+
+        // Inline accesor functions
+        ScalarT& V()
+        {
+            return bus_->V();
+        }
+
+        const ScalarT& V() const
+        {
+            return bus_->V();
+        }
+
+        ScalarT& theta()
+        {
+            return bus_->theta();
+        }
+
+        const ScalarT& theta() const
+        {
+            return bus_->theta();
+        }
+        ScalarT& P()
+        {
+            return bus_->P();
+        }
+
+        const ScalarT& P() const
+        {
+            return bus_->P();
+        }
+
+        ScalarT& Q()
+        {
+            return bus_->Q();
+        }
+
+        const ScalarT& Q() const
+        {
+            return bus_->Q();
+        }
+
+
+    private:
+        real_type omega0_;///< Nominal frequency (assume 2pi60 rad/s)
+        real_type H_;     ///< Inertia constant [s]
+        real_type D_;     ///< Damping constant [pu]
+        real_type Rs_;    ///< stator resistance [pu]
+        real_type Xl_;    ///< stator leakage reactance [pu]
+        real_type Xd_;    ///< direct axis synchronous reactance [pu]
+        real_type Xdp_;   ///< direct axis transient reactance [pu]
+        real_type Xdpp_;  ///< direct axis sub-transient reactance [pu]
+        real_type Xq_;    ///< quadrature axis synchronous reactance [pu]
+        real_type Xqp_;   ///< quadrature axis transient reactance [pu]
+        real_type Xqpp_;  ///< quadrature axis sub-transient reactance [pu]
+        real_type Td0p_;  ///< open circuit direct axis transient time constant [s]
+        real_type Td0pp_; ///< open circuit direct axis sub-transient time constant [s]
+        real_type Tq0p_;  ///< open circuit quadrature axis transient time constant [s]
+        real_type Tq0pp_; ///< open circuit quadrature axis sub-transient time constant [s]
+        real_type S1_;    ///< saturation at 1.0 pu flux
+        real_type S12_;   ///< saturation at 1.2 pu flux
+        real_type X_;
+        real_type Ef_;
+        real_type omega_s_;
+
+        real_type R_;
+        real_type T1_;
+        real_type T2_;
+        real_type T3_;
+        real_type Vmax_;
+        real_type Vmin_;
+        real_type Dt_;
+        real_type offsetGov_;
+
+        ScalarT P0_;
+        ScalarT Q0_;
+
+        bus_type* bus_;
+    };
+
+} // namespace ModelLib
+
+
+#endif // _GOVERNOR_H
\ No newline at end of file
diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/CMakeLists.txt b/ComponentLib/DynamicPhasor/SynchronousMachine/CMakeLists.txt
new file mode 100644
index 0000000..279e0eb
--- /dev/null
+++ b/ComponentLib/DynamicPhasor/SynchronousMachine/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(GENROUwS)
\ No newline at end of file
diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/CMakeLists.txt b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/CMakeLists.txt
new file mode 100644
index 0000000..8fd87d9
--- /dev/null
+++ b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/CMakeLists.txt
@@ -0,0 +1,5 @@
+gridkit_add_library(genrou
+  SOURCES
+    GENROU.cpp
+  OUTPUT_NAME
+    gridkit_genrou)
\ No newline at end of file
diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp
new file mode 100644
index 0000000..e18018e
--- /dev/null
+++ b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.cpp
@@ -0,0 +1,354 @@
+// Made by Paul Moon 6/7/2024
+
+#include <iostream>
+#include <cmath>
+#include <complex>
+
+#include <ComponentLib/PowerFlow/Bus/BaseBus.hpp>
+#include "GENROU.hpp"
+
+namespace ModelLib
+{
+
+    /*!
+     * @brief Constructor for a simple generator model
+     *
+     * Arguments passed to ModelEvaluatorImpl:
+     * - Number of equations = 6 differential + 9 algebraic = 15
+     * - Quadratures and optimization parameters are not being used.
+     */
+    template <class ScalarT, typename IdxT>
+    GENROU<ScalarT, IdxT>::GENROU(bus_type *bus, ScalarT P0, ScalarT Q0)
+        : ModelEvaluatorImpl<ScalarT, IdxT>(15, 0, 0),
+          H_(3.0),                    // Intertia constant
+          D_(0.0),                    // Damping coefficient
+          Xq_(0.5),                   // Machine reactance parameter
+          Xd_(2.1),                   // Machine reactance parameter
+          Xqp_(0.5),                  // Machine reactance parameter
+          Xdp_(0.2),                  // Machine reactance parameter
+          Xqpp_(0.18),                // Machine reactance parameter
+          Xdpp_(0.18),                // Machine reactance parameter
+          Rs_(0.0),                   // Winding resistance
+          Tq0p_(0.75),                // Time constant
+          Td0p_(7.0),                 // Time constant
+          Tq0pp_(0.05),               // Time constant
+          Td0pp_(0.035),              // Time constant
+          Ef_(1.8336),                // Field winding voltage from the excitation system. An exciter would set this value
+          Pm_(1.0),                   // Mechanical power from the prime mover. A governor would set this value
+          X_(0.22),                   // Machine reactance parameter
+          Xl_(0.15),                  // Machine reactance parameter
+          omega_s_(0.0),              // Relative speed deviation. 0 is synchronous, -1 is 0 Hz.
+          omega0_(2.0 * 60.0 * M_PI), // Nominal frequency
+          S1_(0.05),                  // Saturation parameter
+          S12_(0.2),                  // Saturation parameter
+          P0_(1.0),                   // Real power
+          Q0_(0.0),                   // Reactive power
+          bus_(bus)
+    {
+    }
+
+    template <class ScalarT, typename IdxT>
+    GENROU<ScalarT, IdxT>::~GENROU()
+    {
+    }
+
+    /*!
+     * @brief This function will be used to allocate sparse Jacobian matrices.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int GENROU<ScalarT, IdxT>::allocate()
+    {
+        tag_.resize(size_);
+
+        return 0;
+    }
+
+    /**
+     * @brief Initialization of the generator model
+     *
+     * Initialization equations are derived from Adam
+     * Birchfield's paper on power electronics.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int GENROU<ScalarT, IdxT>::initialize()
+    {
+        // Compute values for saturation constants and real/imaginary voltage values
+        double SB = ((1 / 0.2) * (1 / 0.2)) * ((sqrt(S1_) - sqrt(S12_)) * (sqrt(S1_) - sqrt(S12_)));
+        double SA = 1 - sqrt(S1_ / SB);
+        const double Vr = V() * cos(theta());
+        const double Vi = V() * sin(theta());
+
+        const double Ir = P0_ / V();
+        const double Ii = 0;
+
+        // Compute initial guess for the generator voltage phase
+        const double delta = atan((Xq_ * Ir + Vi) / (-Xq_ * Ii + Vr));
+
+        // Compute initial guesses for generator currents and potentials in d-q frame
+        const double Id = Ir * sin(delta) - Ii * cos(delta);
+        const double Iq = Ir * cos(delta) + Ii * sin(delta);
+
+        const double Vq = Vr * cos(delta) + Vi * sin(delta) + Id * Xqpp_ + Iq * Rs_;
+        const double Vd = Vr * sin(delta) - Vi * cos(delta) + Id * Rs_ - Iq * Xqpp_;
+
+        // Compute initial guesses for flux
+        const double Psiqpp = -Vd / (1 + omega_s_);
+        const double Psidpp = Vq / (1 + omega_s_);
+
+        const double Psipp = sqrt((Psiqpp * Psiqpp) + (Psidpp * Psidpp));
+        const double ksat = SB * ((Psipp - SA) * (Psipp - SA));
+
+        // Compute initial guess for Telec. Note that Pmech would normally be changed by a governor.
+        const double Telec = (Psidpp - Id * Xdpp_) * Iq - (Psiqpp - Iq * Xdpp_) * Id;
+        const double Pmech = Telec;
+
+        // Compute remaining initial guesses.
+        const double Psiqp = (1.0 / (1.0 + (Xqp_ - Xqpp_) / (Xqpp_ - Xl_))) * (-Psiqpp * (Xqp_ - Xl_) / (Xqpp_ - Xl_) + Iq * (Xqp_ - Xl_));
+
+        const double Edp = Psiqp - Iq * (Xqp_ - Xl_);
+
+        const double Psidp = -(1.0 / (1.0 + (Xdp_ - Xdpp_) / (Xdpp_ - Xl_))) * (Id * (Xdp_ - Xl_) - Psidpp * (Xdp_ - Xl_) / (Xdpp_ - Xl_));
+
+        const double Eqp = Psidp + Id * (Xdp_ - Xl_);
+
+        y_[0] = delta;
+        y_[1] = omega_s_;
+        y_[2] = Eqp;
+        y_[3] = Psidp;
+        y_[4] = Psiqp;
+        y_[5] = Edp;
+        y_[6] = Psiqpp;
+        y_[7] = Psidpp;
+        y_[8] = Psipp;
+        y_[9] = ksat;
+        y_[10] = Vd;
+        y_[11] = Vq;
+        y_[12] = Telec;
+        y_[13] = Id;
+        y_[14] = Iq;
+
+        yp_[0] = 0.0;
+        yp_[1] = 0.0;
+        yp_[2] = 0.0;
+        yp_[3] = 0.0;
+        yp_[4] = 0.0;
+        yp_[5] = 0.0;
+        yp_[6] = 0.0;
+        yp_[7] = 0.0;
+        yp_[8] = 0.0;
+        yp_[9] = 0.0;
+        yp_[10] = 0.0;
+        yp_[11] = 0.0;
+        yp_[12] = 0.0;
+        yp_[13] = 0.0;
+        yp_[14] = 0.0;
+
+        return 0;
+    }
+
+    /**
+     * \brief Identify differential variables.
+     */
+    template <class ScalarT, typename IdxT>
+    int GENROU<ScalarT, IdxT>::tagDifferentiable()
+    {
+        tag_[0] = true;
+        tag_[1] = true;
+        tag_[2] = true;
+        tag_[3] = true;
+        tag_[4] = true;
+        tag_[5] = true;
+        for (IdxT i = 6; i < size_; ++i)
+        {
+            tag_[i] = false;
+        }
+
+        return 0;
+    }
+
+    /**
+     * @brief Computes residual vector for the generator model.
+     *
+     * Residual equations are given by Adam Birchfield's paper
+     * on power electronics.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int GENROU<ScalarT, IdxT>::evaluateResidual()
+    {
+        // These constants are the same as in the initialization state. Maybe move them to .hpp file?
+        double SB = ((1 / 0.2) * (1 / 0.2)) * ((sqrt(S1_) - sqrt(S12_)) * (sqrt(S1_) - sqrt(S12_)));
+        double SA = 1 - sqrt(S1_ / SB);
+        const double Vr = V() * cos(theta());
+        const double Vi = V() * sin(theta());
+
+        // The first 6 functions are differential equations (involving yp_)
+        f_[0] = -yp_[0] + omega0_ * y_[1];
+        f_[1] = -yp_[1] + (1 / (2 * H_)) * ((Pm_ - D_ * y_[1]) / (1 + y_[1]) - y_[12]);
+        f_[2] = -yp_[2] * Td0p_ + Ef_ - (y_[2] + (Xd_ - Xdp_) * (y_[13] + ((Xdp_ - Xdpp_) / ((Xdp_ - Xl_) * (Xdp_ - Xl_))) * (y_[2] - y_[3] - (Xdp_ - Xl_) * y_[13])) + y_[7] * y_[9]);
+        f_[3] = -yp_[3] + (1 / Td0pp_) * (y_[2] - y_[3] - (Xdp_ - Xl_) * y_[13]);
+        f_[4] = -yp_[4] + (1 / Tq0pp_) * (y_[5] - y_[4] + (Xqp_ - Xl_) * y_[14]);
+        f_[5] = -yp_[5] + (1 / Tq0p_) * (-y_[5] + (((Xq_ - Xl_) / (Xd_ - Xl_)) * y_[6] * y_[9] + (Xq_ - Xqp_) * (y_[14] - ((Xqp_ - Xqpp_) / ((Xqp_ - Xl_) * (Xqp_ - Xl_)) * y_[5] + y_[14]) * (Xqp_ - Xl_) - y_[4])));
+
+        // Remaining equations are all algebraic
+        f_[6] = -y_[6] - y_[4] * ((Xqp_ - Xqpp_) / (Xqp_ - Xl_)) - y_[5] * ((Xqpp_ - Xl_) / (Xqp_ - Xl_));
+        f_[7] = -y_[7] + y_[3] * ((Xdp_ - Xdpp_) / (Xdp_ - Xl_)) + y_[2] * ((Xdpp_ - Xl_) / (Xdp_ - Xl_));
+        f_[8] = -y_[8] + sqrt((y_[7] * y_[7]) + (y_[6] * y_[6]));
+        f_[9] = -y_[9] + SB * ((y_[8] - SA) * (y_[8] - SA));
+        f_[10] = -y_[10] - y_[6] * (1 + y_[1]);
+        f_[11] = -y_[11] + y_[7] * (1 + y_[1]);
+        f_[12] = -y_[12] + y_[14] * (y_[7] - y_[13] * Xdpp_) - y_[13] * (y_[6] - y_[14] * Xdpp_);
+        f_[13] = Vr * sin(y_[0]) - Vi * cos(y_[0]) + y_[13] * Rs_ - y_[14] * Xqpp_;
+        f_[14] = Vr * cos(y_[0]) - Vi * sin(y_[0]) + y_[13] * Xqpp_ - y_[14] * Rs_;
+
+        return 0;
+    }
+
+    /**
+     * @brief Creates the Jacobian using J and Jp.
+     *
+     * Values were calculated by hand and verified with MATLAB.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int GENROU<ScalarT, IdxT>::evaluateJacobian()
+    {
+        J_.zeroMatrix();
+        double SB = ((1 / 0.2) * (1 / 0.2)) * ((sqrt(S1_) - sqrt(S12_)) * (sqrt(S1_) - sqrt(S12_)));
+        double SA = 1 - sqrt(S1_ / SB);
+        const double Vr = V() * cos(theta());
+        const double Vi = V() * sin(theta());
+
+        // These variables are for ease of reading
+        double df2_domega = Pm_ / ((1 + (y_[1]) * (y_[1])) * 2 * H_) - D_ / ((1 + y_[1]) * 2 * H_) + D_ * y_[1] / ((1 + (y_[1]) * (y_[1])) * 2 * H_);
+        double df2_dy12 = -1 / (2 * H_);
+
+        double df3_dId = -(Xd_ - Xdp_) * (1 - (Xdp_ - Xl_) * ((Xdp_ - Xdpp_) / ((Xdp_ - Xl_) * (Xdp_ - Xl_))));
+        double df3_dksat = -y_[7];
+        double df3_dpsippd = -y_[9];
+        double df3_dpsipd = (Xd_ - Xdp_) * (Xdp_ - Xdpp_) / ((Xdp_ - Xl_) * (Xdp_ - Xl_));
+        double df3_dEpq = -1 - (Xd_ - Xdp_) * ((Xdp_ - Xdpp_) / ((Xdp_ - Xl_) * (Xdp_ - Xl_)));
+
+        double df6_dpsippq = y_[9] * ((Xq_ - Xl_) / (Xd_ - Xl_)) / Tq0p_;
+        double df6_dEpd = -(1 / Tq0p_) * (1 + (Xq_ - Xqp_) * ((Xqp_ - Xqpp_) / ((Xqp_ - Xl_) * (Xqp_ - Xl_))));
+        double df6_dpsipq = (Xq_ - Xqp_) * ((Xqp_ - Xqpp_) / ((Xqp_ - Xl_) * (Xqp_ - Xl_))) / Tq0p_;
+        double df6_dksat = y_[6] * ((Xq_ - Xl_) / (Xd_ - Xl_)) / Tq0p_;
+        double df6_dIq = (Xq_ - Xqp_) / Tq0p_ - (Xqp_ - Xl_) * (Xq_ - Xqp_) * ((Xqp_ - Xqpp_) / ((Xqp_ - Xl_) * (Xqp_ - Xl_))) / Tq0p_;
+
+        double df7_dpsipq = -((Xqp_ - Xqpp_) / (Xqp_ - Xl_));
+        double df7_dEpd = -((Xqpp_ - Xl_) / (Xqp_ - Xl_));
+
+        double df8_dpsipd = (Xdp_ - Xdpp_) / (Xdp_ - Xl_);
+        double df8_dEpq = (Xdpp_ - Xl_) / (Xdp_ - Xl_);
+
+        double dTelec_dId = -Xdpp_ * y_[14] - y_[6] + y_[14] * Xdpp_;
+        double dTelec_dIq = y_[7] - y_[13] * Xdpp_ + Xdpp_ * y_[13];
+
+        // Create dF/dy
+        std::vector<IdxT> rcord{0,
+                                1, 1,
+                                2, 2, 2, 2, 2,
+                                3, 3, 3,
+                                4, 4, 4,
+                                5, 5, 5, 5, 5,
+                                6, 6, 6,
+                                7, 7, 7,
+                                8, 8, 8,
+                                9, 9,
+                                10, 10, 10,
+                                11, 11, 11,
+                                12, 12, 12, 12, 12,
+                                13, 13, 13,
+                                14, 14, 14};
+
+        std::vector<IdxT> ccord{1,
+                                1, 12,
+                                2, 3, 7, 9, 13,
+                                2, 3, 13,
+                                4, 5, 14,
+                                4, 5, 6, 9, 14,
+                                4, 5, 6,
+                                2, 3, 7,
+                                6, 7, 8,
+                                8, 9,
+                                1, 6, 10,
+                                1, 7, 11,
+                                6, 7, 12, 13, 14,
+                                0, 13, 14,
+                                0, 13, 14};
+
+        std::vector<ScalarT> vals{omega0_,
+                                  -df2_domega, df2_dy12,
+                                  df3_dEpq, df3_dpsipd, df3_dpsippd, df3_dksat, df3_dId,
+                                  1 / Td0pp_, -1 / Td0pp_, (Xl_ - Xdp_) / Td0pp_,
+                                  -1 / Tq0pp_, 1 / Tq0pp_, (Xqp_ - Xl_) / Tq0pp_,
+                                  df6_dpsipq, df6_dEpd, df6_dpsippq, df6_dksat, df6_dIq,
+                                  df7_dpsipq, df7_dEpd, -1,
+                                  df8_dEpq, df8_dpsipd, -1,
+                                  y_[6] / sqrt((y_[7] * y_[7]) + (y_[6] * y_[6])), y_[7] / sqrt((y_[7] * y_[7]) + (y_[6] * y_[6])), -1,
+                                  2.0 * SB * (y_[8] - SA), -1,
+                                  -y_[6], -(1 + y_[1]), -1,
+                                  y_[7], 1 + y_[1], -1,
+                                  -y_[13], y_[14], -1, dTelec_dId, dTelec_dIq,
+                                  Vr * cos(y_[0]) + Vi * sin(y_[0]), Rs_, -Xqpp_,
+                                  -Vr * sin(y_[0]) + Vi * cos(y_[0]), Xqpp_, Rs_};
+        J_.setValues(rcord, ccord, vals);
+
+        // Create dF/dy'
+        std::vector<IdxT> rcordder{0, 1, 2, 3, 4, 5};
+        std::vector<IdxT> ccordder{0, 1, 2, 3, 4, 5};
+        std::vector<ScalarT> valsder{-1, -1, -Td0p_, -1, -1, -1};
+        COO_Matrix<ScalarT, IdxT> Jacder = COO_Matrix<ScalarT, IdxT>(rcordder, ccordder, valsder, 15, 15);
+
+        // Perform dF/dy + \alpha dF/dy'
+        J_.axpy(alpha_, Jacder);
+        return 0;
+    }
+
+    /**
+     * @brief This is not used, but is implemented to keep SystemModel.hpp happy.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int GENROU<ScalarT, IdxT>::evaluateIntegrand()
+    {
+        return 0;
+    }
+
+    /**
+     * @brief This is not used, but is implemented to keep SystemModel.hpp happy.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int GENROU<ScalarT, IdxT>::initializeAdjoint()
+    {
+        return 0;
+    }
+
+    /**
+     * @brief This is not used, but is implemented to keep SystemModel.hpp happy.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int GENROU<ScalarT, IdxT>::evaluateAdjointResidual()
+    {
+        return 0;
+    }
+
+    /**
+     * @brief This is not used, but is implemented to keep SystemModel.hpp happy.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int GENROU<ScalarT, IdxT>::evaluateAdjointIntegrand()
+    {
+        return 0;
+    }
+
+    // Available template instantiations
+    template class GENROU<double, long int>;
+    template class GENROU<double, size_t>;
+
+} // namespace ModelLib
\ No newline at end of file
diff --git a/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp
new file mode 100644
index 0000000..16ac561
--- /dev/null
+++ b/ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp
@@ -0,0 +1,118 @@
+// Made by Paul Moon 6/7/2024
+
+#ifndef _GENROU_H_
+#define _GENROU_H_
+
+#include <ModelEvaluatorImpl.hpp>
+
+namespace ModelLib
+{
+    template <class ScalarT, typename IdxT> class BaseBus;
+}
+
+namespace ModelLib
+{
+    /*!
+     * @brief Implementation of a GENROU synchronous machine model.
+     *
+     */
+    template  <class ScalarT, typename IdxT>
+    class GENROU : public ModelEvaluatorImpl<ScalarT, IdxT>
+    {
+        using ModelEvaluatorImpl<ScalarT, IdxT>::size_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::nnz_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::time_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::alpha_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::y_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::yp_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::tag_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::f_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::g_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::yB_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::ypB_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::fB_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::gB_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::param_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::param_up_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::param_lo_;
+        using ModelEvaluatorImpl<ScalarT, IdxT>::J_;
+
+        typedef typename ModelEvaluatorImpl<ScalarT, IdxT>::real_type real_type;
+        typedef BaseBus<ScalarT, IdxT> bus_type;
+
+    public:
+        GENROU(BaseBus<ScalarT, IdxT>* bus, ScalarT P0 = 1.0, ScalarT Q0 = 0.0);
+        virtual ~GENROU();
+
+        int allocate();
+        int initialize();
+        int tagDifferentiable();
+        int evaluateResidual();
+        int evaluateJacobian();
+
+        int evaluateIntegrand();
+        int initializeAdjoint();
+        int evaluateAdjointResidual();
+        int evaluateAdjointIntegrand();
+
+        void updateTime(real_type t, real_type a)
+        {
+            time_ = t;
+            alpha_ = a;
+        }
+
+        // Inline accesor functions
+        ScalarT& V()
+        {
+            return bus_->V();
+        }
+
+        const ScalarT& V() const
+        {
+            return bus_->V();
+        }
+
+        ScalarT& theta()
+        {
+            return bus_->theta();
+        }
+
+        const ScalarT& theta() const
+        {
+            return bus_->theta();
+        }
+
+    private:
+        real_type omega0_;///< Nominal frequency (assume 2pi60 rad/s)
+        real_type H_;     ///< Inertia constant [s]
+        real_type D_;     ///< Damping constant [pu]
+        real_type Rs_;    ///< stator resistance [pu]
+        real_type Xl_;    ///< stator leakage reactance [pu]
+        real_type Xd_;    ///< direct axis synchronous reactance [pu]
+        real_type Xdp_;   ///< direct axis transient reactance [pu]
+        real_type Xdpp_;  ///< direct axis sub-transient reactance [pu]
+        real_type Xq_;    ///< quadrature axis synchronous reactance [pu]
+        real_type Xqp_;   ///< quadrature axis transient reactance [pu]
+        real_type Xqpp_;  ///< quadrature axis sub-transient reactance [pu]
+        real_type Td0p_;  ///< open circuit direct axis transient time constant [s]
+        real_type Td0pp_; ///< open circuit direct axis sub-transient time constant [s]
+        real_type Tq0p_;  ///< open circuit quadrature axis transient time constant [s]
+        real_type Tq0pp_; ///< open circuit quadrature axis sub-transient time constant [s]
+        real_type S1_;    ///< saturation at 1.0 pu flux
+        real_type S12_;   ///< saturation at 1.2 pu flux
+        real_type X_;
+
+        real_type Ef_;
+        real_type Pm_;
+        real_type omega_s_;
+
+        ScalarT P0_;
+        ScalarT Q0_;
+
+        bus_type* bus_;
+    };
+
+} // namespace ModelLib
+
+
+#endif // _GENROU_H_
diff --git a/ComponentLib/PowerFlow/Generator4Governor/Generator4Governor.cpp b/ComponentLib/PowerFlow/Generator4Governor/Generator4Governor.cpp
index db26213..11dba5c 100644
--- a/ComponentLib/PowerFlow/Generator4Governor/Generator4Governor.cpp
+++ b/ComponentLib/PowerFlow/Generator4Governor/Generator4Governor.cpp
@@ -157,6 +157,7 @@ int Generator4Governor<ScalarT, IdxT>::initialize()
 
     // Compute generator voltage phase
     const ScalarT delta = atan((Xq_*P0_ - Rs_*Q0_) / (V()*V() + Rs_*P0_ + Xq_*Q0_)) + theta();
+    std::cout << delta << std::endl;
 
     // Compute generator current phase
     const ScalarT phi   = delta - theta() + atan(Q0_/P0_);
@@ -169,7 +170,7 @@ int Generator4Governor<ScalarT, IdxT>::initialize()
 
     // Initialize generator
     y_[offsetGen_ + 0] =  delta;
-    y_[offsetGen_ + 1] =  omega_s_ + 0.2; // <~ this is hack to perturb omega
+    y_[offsetGen_ + 1] =  omega_s_; // <~ this is hack to perturb omega
     y_[offsetGen_ + 2] =  Edp;
     y_[offsetGen_ + 3] =  Eqp;
     y_[offsetGen_ + 4] =  Id;
diff --git a/Examples/CMakeLists.txt b/Examples/CMakeLists.txt
index 67e0da1..2283111 100644
--- a/Examples/CMakeLists.txt
+++ b/Examples/CMakeLists.txt
@@ -73,5 +73,6 @@ if(TARGET SUNDIALS::idas)
     add_subdirectory(GenConstLoad)
     add_subdirectory(GenInfiniteBus)
     add_subdirectory(ParameterEstimation)
+    add_subdirectory(GenrouTest)
   endif()
 endif()
diff --git a/Examples/GenConstLoad/GenConstLoad.cpp b/Examples/GenConstLoad/GenConstLoad.cpp
index fbbc271..f5486a2 100644
--- a/Examples/GenConstLoad/GenConstLoad.cpp
+++ b/Examples/GenConstLoad/GenConstLoad.cpp
@@ -60,6 +60,8 @@
 
 #include <iostream>
 #include <iomanip>
+#include <fstream>
+#include <filesystem>
 
 #include <ComponentLib/PowerFlow/Bus/BusPQ.hpp>
 #include <ComponentLib/PowerFlow/Load/Load.hpp>
@@ -73,6 +75,20 @@
 #include <Solver/Optimization/DynamicConstraint.hpp>
 #include <Utilities/Testing.hpp>
 
+void printOutput(realtype time, const N_Vector yy)
+{
+    static std::ofstream outFile("output.txt");
+
+    realtype* yval = N_VGetArrayPointer_Serial(yy);
+
+    outFile << time;
+    for (size_t j = 0; j < N_VGetLength_Serial(yy); j++)
+    {
+        outFile << " " << yval[j];
+    }
+    outFile << "\n";
+}
+
 int main()
 {
     using namespace ModelLib;
@@ -98,8 +114,45 @@ int main()
     // Create numerical integrator and configure it for the generator model
     Ida<double, size_t>* idas = new Ida<double, size_t>(model);
 
+    // allocate model components
     model->allocate();
 
+    std::cout << "Size: " << model->y().size() << std::endl;
+
+    model->initialize();
+    model->evaluateResidual();
+
+    std::cout << "Verify Initial Residual is Zero: {";
+    for (double i : model->getResidual())
+    {
+        std::cout << i << ", ";
+    }
+    std::cout << "}" << std::endl;
+
+
+    //model->updateTime(0.0, 1.0);
+    //model->evaluateJacobian();
+    //std::cout << "Initial Jacobian with alpha = 1:\n";
+    //model->getJacobian().printMatrix();
+    
+
+    // Create numerical integrator and configure it for the generator model
+    //AnalysisManager::Sundials::Ida<double, size_t>* idas = new AnalysisManager::Sundials::Ida<double, size_t>(model);
+
+    double t_init  = 0.0;
+    double t_final = 5.0;
+    double t_timestep = 0.001;
+
+    idas->setOutputCallback(printOutput);
+
+    idas->configureSimulation();
+    idas->getDefaultInitialCondition();
+    idas->initializeSimulation(t_init, true);
+    idas->runSimulation(t_final, 5000);
+
+
+    /*model->allocate();
+
     double t_init  = 0.0;
     double t_final = 15.0;
 
@@ -171,7 +224,7 @@ int main()
     if(retval != 0)
     {
         std::cout << "The two results differ beyond solver tolerance!\n";
-    }
+    }*/
 
 
     delete idas;
diff --git a/Examples/GenrouTest/CMakeLists.txt b/Examples/GenrouTest/CMakeLists.txt
new file mode 100644
index 0000000..c10f723
--- /dev/null
+++ b/Examples/GenrouTest/CMakeLists.txt
@@ -0,0 +1,9 @@
+add_executable(genroutest GenrouTest.cpp)
+add_executable(genrougovtest GenrouGovTest.cpp)
+target_link_libraries(genroutest GRIDKIT::bus GRIDKIT::genrou GRIDKIT::solvers_opt GRIDKIT::solvers_dyn)
+target_link_libraries(genrougovtest GRIDKIT::bus GRIDKIT::genrou GRIDKIT::solvers_opt GRIDKIT::solvers_dyn GRIDKIT::governor GRIDKIT::load)
+install(TARGETS genroutest DESTINATION bin)
+install(TARGETS genrougovtest DESTINATION bin)
+
+add_test(NAME GenrouTest COMMAND $<TARGET_FILE:genroutest>)
+add_test(NAME GenrouGovTest COMMAND $<TARGET_FILE:genrougovtest>)
\ No newline at end of file
diff --git a/Examples/GenrouTest/GenrouGovTest.cpp b/Examples/GenrouTest/GenrouGovTest.cpp
new file mode 100644
index 0000000..df723e4
--- /dev/null
+++ b/Examples/GenrouTest/GenrouGovTest.cpp
@@ -0,0 +1,102 @@
+// Made by Paul Moon 7/11/2024
+
+
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <filesystem>
+
+#include <ComponentLib/PowerFlow/Bus/BusPQ.hpp>
+#include <ComponentLib/PowerFlow/Load/Load.hpp>
+#include <ComponentLib/DynamicPhasor/Governor/Governor.hpp>
+#include <SystemModel.hpp>
+#include <Solver/Dynamic/Ida.hpp>
+
+#include <Solver/Optimization/DynamicObjective.hpp>
+#include <Solver/Optimization/DynamicConstraint.hpp>
+#include <IpIpoptApplication.hpp>
+#include <IpSolveStatistics.hpp>
+#include <Utilities/Testing.hpp>
+
+void printOutput(realtype time, const N_Vector yy)
+{
+    static std::ofstream outFile("output.txt");
+
+    realtype* yval = N_VGetArrayPointer_Serial(yy);
+
+    outFile << time;
+    for (size_t j = 0; j < N_VGetLength_Serial(yy); j++)
+    {
+        outFile << " " << yval[j];
+    }
+    outFile << "\n";
+}
+
+int main()
+{
+    using namespace ModelLib;
+    using namespace AnalysisManager::Sundials;
+    using namespace AnalysisManager;
+    using namespace GridKit::Testing;
+    
+
+    // Create a bus
+    BaseBus<double, size_t>* bus = new BusPQ<double, size_t>(1.0, 0.0);
+
+    // In GenrouGovTest.cpp
+    ModelEvaluatorImpl<double, size_t>* gen = new Governor<double, size_t>(bus, 1.0, 0.0);
+
+    // Attach load to the bus
+    ModelEvaluatorImpl<double, size_t>* load = new Load<double, size_t>(bus, 1.0, 0.0);
+
+    // Create a system model
+    SystemModel<double, size_t>* model = new SystemModel<double, size_t>();
+    model->addBus(bus);
+    model->addComponent(gen);
+    model->addComponent(load);
+
+    Ida<double, size_t>* idas = new Ida<double, size_t>(model);
+
+    // allocate model components
+    model->allocate();
+
+    std::cout << "Size: " << model->y().size() << std::endl;
+
+    model->initialize();
+    model->evaluateResidual();
+
+    std::cout << "Verify Initial Residual is Zero: {";
+    for (double i : model->getResidual())
+    {
+        std::cout << i << ", ";
+    }
+    std::cout << "}" << std::endl;
+
+
+    //model->updateTime(0.0, 1.0);
+    //model->evaluateJacobian();
+    //std::cout << "Initial Jacobian with alpha = 1:\n";
+    //model->getJacobian().printMatrix();
+    
+
+    // Create numerical integrator and configure it for the generator model
+    //AnalysisManager::Sundials::Ida<double, size_t>* idas = new AnalysisManager::Sundials::Ida<double, size_t>(model);
+
+    double t_init  = 0.0;
+    double t_final = 1.0;
+    double t_timestep = 0.001;
+
+    idas->setOutputCallback(printOutput);
+
+    idas->configureSimulation();
+    idas->getDefaultInitialCondition();
+    idas->initializeSimulation(t_init, true);
+    idas->runSimulation(t_final, 5000);
+
+    delete idas;
+    delete gen;
+    delete load;
+    delete bus;
+    delete model;
+    return 0;
+}
diff --git a/Examples/GenrouTest/GenrouTest.cpp b/Examples/GenrouTest/GenrouTest.cpp
new file mode 100644
index 0000000..8d410a0
--- /dev/null
+++ b/Examples/GenrouTest/GenrouTest.cpp
@@ -0,0 +1,92 @@
+// Made by Paul Moon 6/7/2024
+
+
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <filesystem>
+
+#include <ComponentLib/PowerFlow/Bus/BusSlack.hpp>
+#include <ComponentLib/DynamicPhasor/SynchronousMachine/GENROUwS/GENROU.hpp>
+#include <SystemModel.hpp>
+#include <Solver/Dynamic/Ida.hpp>
+#include <Solver/Optimization/DynamicObjective.hpp>
+#include <Solver/Optimization/DynamicConstraint.hpp>
+
+#include <IpIpoptApplication.hpp>
+#include <IpSolveStatistics.hpp>
+#include <Utilities/Testing.hpp>
+
+void printOutput(realtype time, const N_Vector yy)
+{
+    static std::ofstream outFile("output.txt");
+
+    realtype* yval = N_VGetArrayPointer_Serial(yy);
+
+    outFile << time;
+    for (size_t j = 0; j < N_VGetLength_Serial(yy); j++)
+    {
+        outFile << " " << yval[j];
+    }
+    outFile << "\n";
+}
+
+int main()
+{
+    using namespace ModelLib;
+    using namespace AnalysisManager::Sundials;
+    using namespace AnalysisManager;
+    using namespace GridKit::Testing;
+    
+
+    // Create an infinite bus
+    BaseBus<double, size_t>* bus = new BusSlack<double, size_t>(1.0, 0.0); 
+
+    // Attach a generator to that bus
+    GENROU<double, size_t>* gen = new GENROU<double, size_t>(bus);
+
+    // Create a system model
+    SystemModel<double, size_t>* model = new SystemModel<double, size_t>();
+    model->addBus(bus);
+    model->addComponent(gen);
+
+    // allocate model components
+    model->allocate();
+
+    std::cout << "Size: " << model->y().size() << std::endl;
+
+    model->initialize();
+    model->evaluateResidual();
+
+    std::cout << "Verify Initial Residual is Zero: {";
+    for (double i : model->getResidual())
+    {
+        std::cout << i << ", ";
+    }
+    std::cout << "}" << std::endl;
+
+
+    model->updateTime(0.0, 1.0);
+    model->evaluateJacobian();
+    //std::cout << "Initial Jacobian with alpha = 1:\n";
+    //model->getJacobian().printMatrix();
+    
+
+    // Create numerical integrator and configure it for the generator model
+    AnalysisManager::Sundials::Ida<double, size_t>* idas = new AnalysisManager::Sundials::Ida<double, size_t>(model);
+
+    double t_init  = 0.0;
+    double t_final = 5.0;
+    double t_timestep = 0.001;
+
+    idas->setOutputCallback(printOutput);
+
+    idas->configureSimulation();
+    idas->getDefaultInitialCondition();
+    idas->initializeSimulation(t_init, true);
+    idas->runSimulation(t_final, 5000);
+
+    delete idas;
+    delete model;
+    return 0;
+}
diff --git a/ModelEvaluatorImpl.hpp b/ModelEvaluatorImpl.hpp
index 0f8969f..28b4c81 100644
--- a/ModelEvaluatorImpl.hpp
+++ b/ModelEvaluatorImpl.hpp
@@ -113,7 +113,7 @@ namespace ModelLib
 
         virtual bool hasJacobian()
         {
-            return false;
+            return !J_.isEmpty();
         }
 
         virtual IdxT size_quad()
diff --git a/Solver/Dynamic/Ida.cpp b/Solver/Dynamic/Ida.cpp
index 637d110..5f0c091 100644
--- a/Solver/Dynamic/Ida.cpp
+++ b/Solver/Dynamic/Ida.cpp
@@ -256,8 +256,8 @@ namespace Sundials
         {
             int initType = IDA_Y_INIT;
 
-            if (tag_)
-                initType = IDA_YA_YDP_INIT;
+            /*if (tag_)
+                initType = IDA_YA_YDP_INIT;*/
 
             retval = IDACalcIC(solver_, initType, 0.1);
             checkOutput(retval, "IDACalcIC");
@@ -284,6 +284,11 @@ namespace Sundials
             checkOutput(retval, "IDASolve");
             //printOutput(tout); 
 
+            if (outputCallback_)
+            {
+                outputCallback_(tret, yy_);
+            }
+
             if (retval == IDA_SUCCESS)
             {
                 ++iout;
@@ -448,7 +453,7 @@ namespace Sundials
         checkOutput(retval, "IDASetUserDataB");
 
         /// \todo Need to set max number of steps based on user input!
-        retval = IDASetMaxNumStepsB(solver_, backwardID_, 2000);
+        retval = IDASetMaxNumStepsB(solver_, backwardID_, 50000);
         checkOutput(retval, "IDASetMaxNumSteps");
 
         // Set up linear solver
@@ -791,6 +796,13 @@ namespace Sundials
         }
     }
 
+    template <class ScalarT, typename IdxT>
+    void Ida<ScalarT, IdxT>::setOutputCallback(std::function<void(realtype, const N_Vector)> callback)
+    {
+        outputCallback_ = callback;
+    }
+
+
     // Compiler will prevent building modules with data type incompatible with realtype
     template class Ida<realtype, long int>;
     template class Ida<realtype, int>;
diff --git a/Solver/Dynamic/Ida.hpp b/Solver/Dynamic/Ida.hpp
index a4cea95..bf12b85 100644
--- a/Solver/Dynamic/Ida.hpp
+++ b/Solver/Dynamic/Ida.hpp
@@ -160,6 +160,7 @@ namespace AnalysisManager
             void printOutput(realtype t);
             void printSpecial(realtype t, N_Vector x);
             void printFinalStats();
+            void setOutputCallback(std::function<void(realtype, const N_Vector)> callback);
 
         private:
             static int Residual(realtype t,
@@ -192,6 +193,7 @@ namespace AnalysisManager
             SUNMatrix JacobianMatB_;
             SUNLinearSolver linearSolver_;
             SUNLinearSolver linearSolverB_;
+            std::function<void(realtype, const N_Vector)> outputCallback_;
 
             real_type t_init_;
             real_type t_final_;
diff --git a/SparseMatrix/COO_Matrix.hpp b/SparseMatrix/COO_Matrix.hpp
index 9036085..c729642 100644
--- a/SparseMatrix/COO_Matrix.hpp
+++ b/SparseMatrix/COO_Matrix.hpp
@@ -69,7 +69,7 @@ class COO_Matrix
     std::tuple<Intdx, Intdx> getDimensions();
 
     void printMatrix();
-
+    bool isEmpty();
     
     static void sortSparseCOO(std::vector<Intdx> &rows, std::vector<Intdx> &columns, std::vector<ScalarT> &vals);
 
@@ -563,7 +563,7 @@ inline void COO_Matrix<ScalarT, Intdx>::printMatrix()
     {
         std::cout << "(" << this->row_indexes_[i]
                   << ", " << this->column_indexes_[i]
-                  << ", " << this->values_[i] << ")\n";
+                  << ", " << this->values_[i] << "),\n";
     }
     std::cout << std::flush;
 }
@@ -726,6 +726,16 @@ inline void COO_Matrix<ScalarT, Intdx>::sortSparseCOO(std::vector<Intdx> &rows,
     }
 }
 
+/**
+ * @brief Return if there are no values, meaning that the matrix is empty.
+ * 
+ */
+template <class ScalarT, typename Intdx>
+inline bool COO_Matrix<ScalarT, Intdx>::isEmpty() 
+{
+    return values_.empty();
+}
+
 template <class ScalarT, typename Intdx>
 inline COO_Matrix<ScalarT, Intdx>::COO_Matrix(std::vector<Intdx> r, std::vector<Intdx> c, std::vector<ScalarT> v, Intdx m, Intdx n)
 {
diff --git a/SystemModel.hpp b/SystemModel.hpp
index 14b7e4e..f2c4587 100644
--- a/SystemModel.hpp
+++ b/SystemModel.hpp
@@ -107,6 +107,7 @@ class SystemModel : public ModelEvaluatorImpl<ScalarT, IdxT>
     using ModelEvaluatorImpl<ScalarT, IdxT>::param_;
     using ModelEvaluatorImpl<ScalarT, IdxT>::param_up_;
     using ModelEvaluatorImpl<ScalarT, IdxT>::param_lo_;
+    using ModelEvaluatorImpl<ScalarT, IdxT>::J_;
 
 public:
     /**
@@ -117,7 +118,7 @@ class SystemModel : public ModelEvaluatorImpl<ScalarT, IdxT>
         // Set system model tolerances
         rtol_ = 1e-7;
         atol_ = 1e-9;
-        this->max_steps_=2000;
+        this->max_steps_=5000;
     }
 
     /**
@@ -192,7 +193,12 @@ class SystemModel : public ModelEvaluatorImpl<ScalarT, IdxT>
      */
     bool hasJacobian() 
     {
-        return false;
+        bool hasJac;
+        for (const auto &component : components_)
+        {
+            hasJac *= component->hasJacobian();
+        }
+        return false; //return hasJac;
     }
 
     /**
@@ -385,12 +391,33 @@ class SystemModel : public ModelEvaluatorImpl<ScalarT, IdxT>
     /**
      * @brief Evaluate system Jacobian.
      *
-     * @todo Need to implement Jacobian. For now, using finite difference
-     * approximation provided by IDA. This works for dense Jacobian matrix
-     * only.
      *
      */
-    int evaluateJacobian(){return 0;}
+    int evaluateJacobian()
+    {
+        IdxT varOffset = 0;
+        IdxT optOffset = 0;
+        for(const auto& component : components_) {
+            for(IdxT j=0; j<component->size(); ++j)
+            {
+                component->y()[j]  = y_[varOffset + j];
+                component->yp()[j] = yp_[varOffset + j];
+            }
+            varOffset += component->size();
+
+            for(IdxT j=0; j<component->size_opt(); ++j)
+            {
+                component->param()[j] = param_[optOffset + j];
+            }
+            optOffset += component->size_opt();
+
+            component->evaluateJacobian();
+        }
+        for(const auto& component : components_) {
+            J_  = component->getJacobian();
+        }
+        return 0;
+    }
 
     /**
      * @brief Evaluate integrands for the system quadratures.