diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp
index 1c8508d1..f70c61c8 100644
--- a/examples/PhasorDynamics/Example1/example1.cpp
+++ b/examples/PhasorDynamics/Example1/example1.cpp
@@ -91,12 +91,10 @@ int main()
   // push it into output, which is updated outside the callback.
   auto output_cb = [&](double t)
   {
-    std::vector<double>& yval = sys.y();
-
     output.push_back({.ti = t,
-                      .Vr = yval[0],
-                      .Vi = yval[1],
-                      .dw = yval[3]});
+                      .Vr = bus1.Vr(),
+                      .Vi = bus1.Vi(),
+                      .dw = gen.y()[1]});
   };
 
   // The above lambda is equivalent to writing
diff --git a/examples/PhasorDynamics/Example2/example2.cpp b/examples/PhasorDynamics/Example2/example2.cpp
index fd580dce..fdab584f 100644
--- a/examples/PhasorDynamics/Example2/example2.cpp
+++ b/examples/PhasorDynamics/Example2/example2.cpp
@@ -111,13 +111,11 @@ int main()
 
   auto output_cb = [&](real_type t)
   {
-    std::vector<double>& yval = sys.y();
-
     output.push_back({.t         = t,
-                      .gen2speed = 1 + yval[5],
-                      .gen3speed = 1 + yval[26],
-                      .v2mag     = sqrt(yval[0] * yval[0] + yval[1] * yval[1]),
-                      .v3mag     = sqrt(yval[2] * yval[2] + yval[3] * yval[3])});
+                      .gen2speed = 1 + gen2.y()[1],
+                      .gen3speed = 1 + gen3.y()[1],
+                      .v2mag     = sqrt(bus2.Vr() * bus2.Vr() + bus2.Vi() * bus2.Vi()),
+                      .v3mag     = sqrt(bus3.Vr() * bus3.Vr() + bus3.Vi() * bus3.Vi())});
   };
 
   /* Set up simulation */
diff --git a/src/Model/Evaluator.hpp b/src/Model/Evaluator.hpp
index b4a447dd..3f94153e 100644
--- a/src/Model/Evaluator.hpp
+++ b/src/Model/Evaluator.hpp
@@ -56,6 +56,12 @@ namespace GridKit
       virtual void setTolerances(real_type& rtol, real_type& atol) const = 0;
       virtual void setMaxSteps(IdxT& msa) const                          = 0;
 
+      /// @brief Let the `Evaluator` know that its state has been updated externally (for example by an integrator),
+      ///        in case it needs to keep other things up to date based on its internal state.
+      virtual void update()
+      {
+      }
+
       virtual std::vector<ScalarT>&       y()       = 0;
       virtual const std::vector<ScalarT>& y() const = 0;
 
diff --git a/src/Model/PhasorDynamics/SystemModel.hpp b/src/Model/PhasorDynamics/SystemModel.hpp
index 03347b22..e2d458b7 100644
--- a/src/Model/PhasorDynamics/SystemModel.hpp
+++ b/src/Model/PhasorDynamics/SystemModel.hpp
@@ -263,42 +263,17 @@ namespace GridKit
        */
       int evaluateResidual()
       {
-        // Update variables
-        IdxT varOffset = 0;
-        IdxT optOffset = 0;
+        // No need to notify children yet, since we're about to
+        // evaluate them
+        updateChildren<false>();
+
         for (const auto& bus : buses_)
         {
-          for (IdxT j = 0; j < bus->size(); ++j)
-          {
-            bus->y()[j]  = y_[varOffset + j];
-            bus->yp()[j] = yp_[varOffset + j];
-          }
-          varOffset += bus->size();
-
-          for (IdxT j = 0; j < bus->sizeParams(); ++j)
-          {
-            bus->param()[j] = param_[optOffset + j];
-          }
-          optOffset += bus->sizeParams();
-
           bus->evaluateResidual();
         }
 
         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->sizeParams(); ++j)
-          {
-            component->param()[j] = param_[optOffset + j];
-          }
-          optOffset += component->sizeParams();
-
           component->evaluateResidual();
         }
 
@@ -343,42 +318,17 @@ namespace GridKit
        */
       int evaluateIntegrand()
       {
-        // Update variables
-        IdxT varOffset = 0;
-        IdxT optOffset = 0;
+        // No need to notify children yet, since we're about to
+        // evaluate them
+        updateChildren<false>();
+
         for (const auto& bus : buses_)
         {
-          for (IdxT j = 0; j < bus->size(); ++j)
-          {
-            bus->y()[j]  = y_[varOffset + j];
-            bus->yp()[j] = yp_[varOffset + j];
-          }
-          varOffset += bus->size();
-
-          for (IdxT j = 0; j < bus->sizeParams(); ++j)
-          {
-            bus->param()[j] = param_[optOffset + j];
-          }
-          optOffset += bus->sizeParams();
-
           bus->evaluateIntegrand();
         }
 
         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->sizeParams(); ++j)
-          {
-            component->param()[j] = param_[optOffset + j];
-          }
-          optOffset += component->sizeParams();
-
           component->evaluateIntegrand();
         }
 
@@ -413,45 +363,11 @@ namespace GridKit
        */
       int initializeAdjoint()
       {
-        IdxT offset    = 0;
-        IdxT optOffset = 0;
-
-        // Update bus variables and optimization parameters
-        for (const auto& bus : buses_)
-        {
-          for (IdxT j = 0; j < bus->size(); ++j)
-          {
-            bus->y()[j]  = y_[offset + j];
-            bus->yp()[j] = yp_[offset + j];
-          }
-          offset += bus->size();
-
-          for (IdxT j = 0; j < bus->sizeParams(); ++j)
-          {
-            bus->param()[j] = param_[optOffset + j];
-          }
-          optOffset += bus->sizeParams();
-        }
-
-        // Update component variables and optimization parameters
-        for (const auto& component : components_)
-        {
-          for (IdxT j = 0; j < component->size(); ++j)
-          {
-            component->y()[j]  = y_[offset + j];
-            component->yp()[j] = yp_[offset + j];
-          }
-          offset += component->size();
-
-          for (IdxT j = 0; j < component->sizeParams(); ++j)
-          {
-            component->param()[j] = param_[optOffset + j];
-          }
-          optOffset += component->sizeParams();
-        }
+        // No need to notify children yet, since we're about to
+        // initialize them
+        updateChildren();
 
-        // Reset counter
-        offset = 0;
+        size_t offset = 0;
 
         // Initialize bus adjoints
         for (const auto& bus : buses_)
@@ -491,45 +407,9 @@ namespace GridKit
        */
       int evaluateAdjointResidual()
       {
-        IdxT varOffset = 0;
-        IdxT optOffset = 0;
-
-        // Update variables in component models
-        for (const auto& bus : buses_)
-        {
-          for (IdxT j = 0; j < bus->size(); ++j)
-          {
-            bus->y()[j]   = y_[varOffset + j];
-            bus->yp()[j]  = yp_[varOffset + j];
-            bus->yB()[j]  = yB_[varOffset + j];
-            bus->ypB()[j] = ypB_[varOffset + j];
-          }
-          varOffset += bus->size();
-
-          for (IdxT j = 0; j < bus->sizeParams(); ++j)
-          {
-            bus->param()[j] = param_[optOffset + j];
-          }
-          optOffset += bus->sizeParams();
-        }
-
-        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];
-            component->yB()[j]  = yB_[varOffset + j];
-            component->ypB()[j] = ypB_[varOffset + j];
-          }
-          varOffset += component->size();
-
-          for (IdxT j = 0; j < component->sizeParams(); ++j)
-          {
-            component->param()[j] = param_[optOffset + j];
-          }
-          optOffset += component->sizeParams();
-        }
+        // No need to notify children yet, since we're about to
+        // evaluate them
+        updateChildrenAdjoint<false>();
 
         for (const auto& bus : buses_)
         {
@@ -575,44 +455,9 @@ namespace GridKit
        */
       int evaluateAdjointIntegrand()
       {
-        // First, update variables
-        IdxT varOffset = 0;
-        IdxT optOffset = 0;
-        for (const auto& bus : buses_)
-        {
-          for (IdxT j = 0; j < bus->size(); ++j)
-          {
-            bus->y()[j]   = y_[varOffset + j];
-            bus->yp()[j]  = yp_[varOffset + j];
-            bus->yB()[j]  = yB_[varOffset + j];
-            bus->ypB()[j] = ypB_[varOffset + j];
-          }
-          varOffset += bus->size();
-
-          for (IdxT j = 0; j < bus->sizeParams(); ++j)
-          {
-            bus->param()[j] = param_[optOffset + j];
-          }
-          optOffset += bus->sizeParams();
-        }
-
-        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];
-            component->yB()[j]  = yB_[varOffset + j];
-            component->ypB()[j] = ypB_[varOffset + j];
-          }
-          varOffset += component->size();
-
-          for (IdxT j = 0; j < component->sizeParams(); ++j)
-          {
-            component->param()[j] = param_[optOffset + j];
-          }
-          optOffset += component->sizeParams();
-        }
+        // No need to notify children yet, since we're about to
+        // evaluate them
+        updateChildrenAdjoint<false>();
 
         // Evaluate integrand and update global vector
         for (const auto& component : components_)
@@ -648,10 +493,113 @@ namespace GridKit
         components_.push_back(component);
       }
 
+      void update() override
+      {
+        updateChildrenAdjoint<true>();
+      }
+
     private:
       std::vector<bus_type*>       buses_;
       std::vector<component_type*> components_;
 
+      /// @brief Whenever the internal state of a system is updated, the system
+      ///        needs to reflect those changes in its "children" aka buses and
+      ///        components. `updateChildren` will only update `y`, `yp`, and
+      ///        `param` - call when you know only the internal state of those
+      ///        variables has been updated.
+      /// @tparam NOTIFY_CHILDREN Whether or not buses and components should
+      ///         be notified through their `update()` method that their state
+      ///         has been updated. If you are planning to notify them otherwise,
+      ///         set this to `false`.
+      template <bool NOTIFY_CHILDREN = true>
+      void updateChildren()
+      {
+        // Update variables
+        IdxT varOffset = 0;
+        IdxT optOffset = 0;
+        for (const auto& bus : buses_)
+        {
+          for (IdxT j = 0; j < bus->size(); ++j)
+          {
+            bus->y()[j]  = y_[varOffset + j];
+            bus->yp()[j] = yp_[varOffset + j];
+          }
+          varOffset += bus->size();
+
+          for (IdxT j = 0; j < bus->sizeParams(); ++j)
+          {
+            bus->param()[j] = param_[optOffset + j];
+          }
+          optOffset += bus->sizeParams();
+
+          if constexpr (NOTIFY_CHILDREN)
+            bus->update();
+        }
+
+        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->sizeParams(); ++j)
+          {
+            component->param()[j] = param_[optOffset + j];
+          }
+          optOffset += component->sizeParams();
+
+          if constexpr (NOTIFY_CHILDREN)
+            component->update();
+        }
+      }
+
+      /// @brief Whenever the internal state of a system is updated, the system
+      ///        needs to reflect those changes in its "children" aka buses and
+      ///        components. `updateChildrenAdjoint` will update all internal
+      ///        variables of this system's children, including adjoint variables.
+      /// @tparam NOTIFY_CHILDREN Whether or not buses and components should
+      ///         be notified through their `update()` method that their state
+      ///         has been updated. If you are planning to notify them otherwise,
+      ///         set this to `false`.
+      template <bool NOTIFY_CHILDREN = true>
+      void updateChildrenAdjoint()
+      {
+        // Update all non-adjoint state.
+        // No need to notify children that they've been updated yet,
+        // because we still have more updating to do.
+        updateChildren<false>();
+
+        IdxT varOffset = 0;
+        IdxT optOffset = 0;
+        for (const auto& bus : buses_)
+        {
+          for (IdxT j = 0; j < bus->size(); ++j)
+          {
+            bus->yB()[j]  = yB_[varOffset + j];
+            bus->ypB()[j] = ypB_[varOffset + j];
+          }
+          varOffset += bus->size();
+
+          if constexpr (NOTIFY_CHILDREN)
+            bus->update();
+        }
+
+        for (const auto& component : components_)
+        {
+          for (IdxT j = 0; j < component->size(); ++j)
+          {
+            component->yB()[j]  = yB_[varOffset + j];
+            component->ypB()[j] = ypB_[varOffset + j];
+          }
+          varOffset += component->size();
+
+          if constexpr (NOTIFY_CHILDREN)
+            component->update();
+        }
+      }
     }; // class SystemModel
 
   } // namespace PhasorDynamics
diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp
index 7cce626d..3d4ba474 100644
--- a/src/Solver/Dynamic/Ida.cpp
+++ b/src/Solver/Dynamic/Ida.cpp
@@ -231,6 +231,7 @@ namespace AnalysisManager
           model_->updateTime(tret, 0.0);
           copyVec(yy_, model_->y());
           copyVec(yp_, model_->yp());
+          model_->update();
 
           (*step_callback)(tret);
         }
@@ -246,6 +247,7 @@ namespace AnalysisManager
       model_->updateTime(tf, 0.0);
       copyVec(yy_, model_->y());
       copyVec(yp_, model_->yp());
+      model_->update();
 
       // std::cout << "\n";
       return retval;