diff --git a/GridKit/LinearAlgebra/Vector/CMakeLists.txt b/GridKit/LinearAlgebra/Vector/CMakeLists.txt
index a32b0cade..5125498ee 100644
--- a/GridKit/LinearAlgebra/Vector/CMakeLists.txt
+++ b/GridKit/LinearAlgebra/Vector/CMakeLists.txt
@@ -1,5 +1,5 @@
 gridkit_add_library(
   dense_vector
   SOURCES Vector.cpp
-  HEADERS Vector.hpp
+  HEADERS Vector.hpp VectorView.hpp
   LINK_LIBRARIES GridKit::utilities_logger)
diff --git a/GridKit/LinearAlgebra/Vector/VectorView.hpp b/GridKit/LinearAlgebra/Vector/VectorView.hpp
new file mode 100644
index 000000000..314a81625
--- /dev/null
+++ b/GridKit/LinearAlgebra/Vector/VectorView.hpp
@@ -0,0 +1,78 @@
+#pragma once
+
+#include <cassert>
+#include <cstddef>
+
+namespace GridKit
+{
+  namespace LinearAlgebra
+  {
+    /**
+     * @brief Non-owning view of contiguous vector data.
+     *
+     * VectorView does not allocate, free, or synchronize memory. It only
+     * provides indexed access to storage owned elsewhere.
+     */
+    template <typename ScalarT, typename IdxT>
+    class VectorView
+    {
+    public:
+      VectorView() = default;
+
+      VectorView(ScalarT* data, std::size_t size)
+        : data_(data),
+          size_(size)
+      {
+        assert(data_ != nullptr || size_ == 0);
+      }
+
+      void bind(ScalarT* data, std::size_t size)
+      {
+        assert(data != nullptr || size == 0);
+        data_ = data;
+        size_ = size;
+      }
+
+      ScalarT* data()
+      {
+        return data_;
+      }
+
+      const ScalarT* data() const
+      {
+        return data_;
+      }
+
+      std::size_t size() const
+      {
+        return size_;
+      }
+
+      IdxT getSize() const
+      {
+        return static_cast<IdxT>(size_);
+      }
+
+      bool empty() const
+      {
+        return size_ == 0;
+      }
+
+      ScalarT& operator[](std::size_t i)
+      {
+        assert(i < size_);
+        return data_[i];
+      }
+
+      const ScalarT& operator[](std::size_t i) const
+      {
+        assert(i < size_);
+        return data_[i];
+      }
+
+    private:
+      ScalarT*    data_{nullptr};
+      std::size_t size_{0};
+    };
+  } // namespace LinearAlgebra
+} // namespace GridKit
diff --git a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp
index a6cd306b2..0abf914e5 100644
--- a/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Branch/BranchEnzyme.cpp
@@ -37,12 +37,12 @@ namespace GridKit
                                      ScalarT,
                                      IdxT>::eval(this,
                                                  static_cast<size_t>(bus1_->size()),
-                                                 static_cast<size_t>((bus1_->y()).size()),
+                                                 static_cast<size_t>(bus1_->size()),
                                                  (bus1_->getResidualIndices()).data(),
                                                  (bus1_->getVariableIndices()).data(),
                                                  y_.data(),
                                                  yp_.data(),
-                                                 (bus1_->y()).data(),
+                                                 bus1_->yData(),
                                                  J_rows_buffer_,
                                                  J_cols_buffer_,
                                                  J_vals_buffer_,
@@ -54,12 +54,12 @@ namespace GridKit
                                      ScalarT,
                                      IdxT>::eval(this,
                                                  static_cast<size_t>(bus2_->size()),
-                                                 static_cast<size_t>((bus2_->y()).size()),
+                                                 static_cast<size_t>(bus2_->size()),
                                                  (bus2_->getResidualIndices()).data(),
                                                  (bus2_->getVariableIndices()).data(),
                                                  y_.data(),
                                                  yp_.data(),
-                                                 (bus2_->y()).data(),
+                                                 bus2_->yData(),
                                                  J_rows_buffer_,
                                                  J_cols_buffer_,
                                                  J_vals_buffer_,
@@ -71,12 +71,12 @@ namespace GridKit
                                      ScalarT,
                                      IdxT>::eval(this,
                                                  static_cast<size_t>(bus1_->size()),
-                                                 static_cast<size_t>((bus2_->y()).size()),
+                                                 static_cast<size_t>(bus2_->size()),
                                                  (bus1_->getResidualIndices()).data(),
                                                  (bus2_->getVariableIndices()).data(),
                                                  y_.data(),
                                                  yp_.data(),
-                                                 (bus2_->y()).data(),
+                                                 bus2_->yData(),
                                                  J_rows_buffer_,
                                                  J_cols_buffer_,
                                                  J_vals_buffer_,
@@ -89,12 +89,12 @@ namespace GridKit
                                      ScalarT,
                                      IdxT>::eval(this,
                                                  static_cast<size_t>(bus2_->size()),
-                                                 static_cast<size_t>((bus1_->y()).size()),
+                                                 static_cast<size_t>(bus1_->size()),
                                                  (bus2_->getResidualIndices()).data(),
                                                  (bus1_->getVariableIndices()).data(),
                                                  y_.data(),
                                                  yp_.data(),
-                                                 (bus1_->y()).data(),
+                                                 bus1_->yData(),
                                                  J_rows_buffer_,
                                                  J_cols_buffer_,
                                                  J_vals_buffer_,
diff --git a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp
index 534af06ed..b3a747085 100644
--- a/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp
@@ -95,9 +95,7 @@ namespace GridKit
       size_t size = static_cast<size_t>(size_);
 
       // Resize component model data
-      f_.resize(size);
-      y_.resize(size);
-      yp_.resize(size);
+      this->allocateState(size);
       tag_.resize(size);
       variable_indices_.resize(size);
       residual_indices_.resize(size);
diff --git a/GridKit/Model/PhasorDynamics/BusBase.hpp b/GridKit/Model/PhasorDynamics/BusBase.hpp
index f6591322c..bfe298e73 100644
--- a/GridKit/Model/PhasorDynamics/BusBase.hpp
+++ b/GridKit/Model/PhasorDynamics/BusBase.hpp
@@ -1,5 +1,7 @@
 #pragma once
 
+#include <cassert>
+#include <cstddef>
 #include <memory>
 #include <set>
 #include <vector>
@@ -7,6 +9,7 @@
 #include <GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp>
 #include <GridKit/CommonMath.hpp>
 #include <GridKit/Constants.hpp>
+#include <GridKit/LinearAlgebra/Vector/VectorView.hpp>
 #include <GridKit/Model/Evaluator.hpp>
 #include <GridKit/Model/PhasorDynamics/Bus/BusData.hpp>
 #include <GridKit/Model/VariableMonitor.hpp>
@@ -27,12 +30,13 @@ namespace GridKit
     class BusBase : public Model::Evaluator<scalar_type, index_type>
     {
     public:
-      using ScalarT  = scalar_type;
-      using IdxT     = index_type;
-      using RealT    = typename Model::Evaluator<ScalarT, IdxT>::RealT;
-      using MatrixT  = typename Model::Evaluator<ScalarT, IdxT>::MatrixT;
-      using BusTypeT = typename BusData<RealT, IdxT>::BusType;
-      using MonitorT = Model::VariableMonitor<BusBase, BusData>;
+      using ScalarT     = scalar_type;
+      using IdxT        = index_type;
+      using RealT       = typename Model::Evaluator<ScalarT, IdxT>::RealT;
+      using MatrixT     = typename Model::Evaluator<ScalarT, IdxT>::MatrixT;
+      using BusTypeT    = typename BusData<RealT, IdxT>::BusType;
+      using MonitorT    = Model::VariableMonitor<BusBase, BusData>;
+      using VectorViewT = LinearAlgebra::VectorView<ScalarT, IdxT>;
 
       BusBase() = default;
 
@@ -66,22 +70,22 @@ namespace GridKit
 
       std::vector<ScalarT>& y() override
       {
-        return y_;
+        return y_storage_;
       }
 
       const std::vector<ScalarT>& y() const override
       {
-        return y_;
+        return y_storage_;
       }
 
       std::vector<ScalarT>& yp() override
       {
-        return yp_;
+        return yp_storage_;
       }
 
       const std::vector<ScalarT>& yp() const override
       {
-        return yp_;
+        return yp_storage_;
       }
 
       std::vector<bool>& tag() override
@@ -96,12 +100,37 @@ namespace GridKit
 
       std::vector<ScalarT>& getResidual() override
       {
-        return f_;
+        return f_storage_;
       }
 
       const std::vector<ScalarT>& getResidual() const override
       {
-        return f_;
+        return f_storage_;
+      }
+
+      ScalarT* yData()
+      {
+        return y_.data();
+      }
+
+      const ScalarT* yData() const
+      {
+        return y_.data();
+      }
+
+      void bindStateViews(VectorViewT y, VectorViewT yp, VectorViewT f)
+      {
+        assert(y.size() == yp.size());
+        assert(y.size() == f.size());
+
+        state_views_bound_ = true;
+        y_storage_.clear();
+        yp_storage_.clear();
+        f_storage_.clear();
+
+        y_  = y;
+        yp_ = yp;
+        f_  = f;
       }
 
       MatrixT& getJacobian() override
@@ -146,6 +175,27 @@ namespace GridKit
         return residual_indices_;
       }
 
+    protected:
+      void allocateState(std::size_t size)
+      {
+        if (state_views_bound_)
+        {
+          assert(y_.size() == size);
+          assert(yp_.size() == size);
+          assert(f_.size() == size);
+          return;
+        }
+
+        f_storage_.resize(size);
+        y_storage_.resize(size);
+        yp_storage_.resize(size);
+
+        f_.bind(size == 0 ? nullptr : f_storage_.data(), size);
+        y_.bind(size == 0 ? nullptr : y_storage_.data(), size);
+        yp_.bind(size == 0 ? nullptr : yp_storage_.data(), size);
+      }
+
+    public:
       /// Pure virtual function, returns bus type (DEFAULT or SLACK).
       virtual BusTypeT BusType() const
       {
@@ -188,12 +238,18 @@ namespace GridKit
       /// Global (system-level) residual indices
       std::vector<IdxT> residual_indices_;
 
-      std::vector<ScalarT> y_;
-      std::vector<ScalarT> yp_;
+      std::vector<ScalarT> y_storage_;
+      std::vector<ScalarT> yp_storage_;
       std::vector<bool>    tag_;
-      std::vector<ScalarT> f_;
+      std::vector<ScalarT> f_storage_;
       std::vector<ScalarT> g_;
 
+      VectorViewT y_;
+      VectorViewT yp_;
+      VectorViewT f_;
+
+      bool state_views_bound_{false};
+
       MatrixT J_;
       IdxT*   J_rows_buffer_{nullptr};
       IdxT*   J_cols_buffer_{nullptr};
diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp
index 9d3c45854..fc71ab139 100644
--- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultEnzyme.cpp
@@ -41,7 +41,7 @@ namespace GridKit
                                                    (bus_->getVariableIndices()).data(),
                                                    y_.data(),
                                                    yp_.data(),
-                                                   (bus_->y()).data(),
+                                                   bus_->yData(),
                                                    J_rows_buffer_,
                                                    J_cols_buffer_,
                                                    J_vals_buffer_,
diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp
index 985c4eae1..e6633fd52 100644
--- a/GridKit/Model/PhasorDynamics/Component.hpp
+++ b/GridKit/Model/PhasorDynamics/Component.hpp
@@ -1,9 +1,12 @@
 #pragma once
 
+#include <cassert>
+#include <cstddef>
 #include <vector>
 
 #include <GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp>
 #include <GridKit/CommonMath.hpp>
+#include <GridKit/LinearAlgebra/Vector/VectorView.hpp>
 #include <GridKit/Model/Evaluator.hpp>
 #include <GridKit/Utilities/Errors.hpp>
 #include <GridKit/Utilities/Logger/Logger.hpp>
@@ -21,8 +24,9 @@ namespace GridKit
     class Component : public Model::Evaluator<ScalarT, IdxT>
     {
     public:
-      using RealT   = typename Model::Evaluator<ScalarT, IdxT>::RealT;
-      using MatrixT = typename Model::Evaluator<ScalarT, IdxT>::MatrixT;
+      using RealT       = typename Model::Evaluator<ScalarT, IdxT>::RealT;
+      using MatrixT     = typename Model::Evaluator<ScalarT, IdxT>::MatrixT;
+      using VectorViewT = LinearAlgebra::VectorView<ScalarT, IdxT>;
 
       Component() = default;
 
@@ -64,22 +68,22 @@ namespace GridKit
 
       std::vector<ScalarT>& y() override
       {
-        return y_;
+        return y_storage_;
       }
 
       const std::vector<ScalarT>& y() const override
       {
-        return y_;
+        return y_storage_;
       }
 
       std::vector<ScalarT>& yp() override
       {
-        return yp_;
+        return yp_storage_;
       }
 
       const std::vector<ScalarT>& yp() const override
       {
-        return yp_;
+        return yp_storage_;
       }
 
       std::vector<bool>& tag() override
@@ -94,12 +98,37 @@ namespace GridKit
 
       std::vector<ScalarT>& getResidual() override
       {
-        return f_;
+        return f_storage_;
       }
 
       const std::vector<ScalarT>& getResidual() const override
       {
-        return f_;
+        return f_storage_;
+      }
+
+      ScalarT* yData()
+      {
+        return y_.data();
+      }
+
+      const ScalarT* yData() const
+      {
+        return y_.data();
+      }
+
+      void bindStateViews(VectorViewT y, VectorViewT yp, VectorViewT f)
+      {
+        assert(y.size() == yp.size());
+        assert(y.size() == f.size());
+
+        state_views_bound_ = true;
+        y_storage_.clear();
+        yp_storage_.clear();
+        f_storage_.clear();
+
+        y_  = y;
+        yp_ = yp;
+        f_  = f;
       }
 
       MatrixT& getJacobian() override
@@ -144,6 +173,27 @@ namespace GridKit
         return residual_indices_;
       }
 
+    protected:
+      void allocateState(std::size_t size)
+      {
+        if (state_views_bound_)
+        {
+          assert(y_.size() == size);
+          assert(yp_.size() == size);
+          assert(f_.size() == size);
+          return;
+        }
+
+        f_storage_.resize(size);
+        y_storage_.resize(size);
+        yp_storage_.resize(size);
+
+        f_.bind(size == 0 ? nullptr : f_storage_.data(), size);
+        y_.bind(size == 0 ? nullptr : y_storage_.data(), size);
+        yp_.bind(size == 0 ? nullptr : yp_storage_.data(), size);
+      }
+
+    public:
       /// @todo Remove this method. It should be part of DynamicSolver class.
       bool hasJacobian() override
       {
@@ -183,12 +233,18 @@ namespace GridKit
       /// Global (system-level) residual indices
       std::vector<IdxT> residual_indices_;
 
-      std::vector<ScalarT> y_;
-      std::vector<ScalarT> yp_;
+      std::vector<ScalarT> y_storage_;
+      std::vector<ScalarT> yp_storage_;
       std::vector<bool>    tag_;
-      std::vector<ScalarT> f_;
+      std::vector<ScalarT> f_storage_;
       std::vector<ScalarT> g_;
 
+      VectorViewT y_;
+      VectorViewT yp_;
+      VectorViewT f_;
+
+      bool state_views_bound_{false};
+
       MatrixT J_;
       IdxT*   J_rows_buffer_{nullptr};
       IdxT*   J_cols_buffer_{nullptr};
diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp
index 31ef04039..80c14fc32 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Enzyme.cpp
@@ -63,7 +63,7 @@ namespace GridKit
                                                    (bus_->getVariableIndices()).data(),
                                                    y_.data(),
                                                    yp_.data(),
-                                                   (bus_->y()).data(),
+                                                   bus_->yData(),
                                                    ws_.data(),
                                                    J_rows_buffer_,
                                                    J_cols_buffer_,
diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp
index c055ad6b4..81b088d5a 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp
@@ -111,9 +111,7 @@ namespace GridKit
       {
         // Resize component model data
         auto size = static_cast<size_t>(size_); // avoid compiler warnings
-        f_.resize(size);
-        y_.resize(size);
-        yp_.resize(size);
+        this->allocateState(size);
         tag_.resize(size);
         variable_indices_.resize(size);
         residual_indices_.resize(size);
diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp
index 2cc959481..40e3c5b4f 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiEnzyme.cpp
@@ -56,7 +56,7 @@ namespace GridKit
                                                    (bus_->getVariableIndices()).data(),
                                                    y_.data(),
                                                    yp_.data(),
-                                                   (bus_->y()).data(),
+                                                   bus_->yData(),
                                                    ws_.data(),
                                                    J_rows_buffer_,
                                                    J_cols_buffer_,
diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp
index e0de70fea..cfb27d956 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp
@@ -58,9 +58,7 @@ namespace GridKit
       int SexsPti<scalar_type, index_type>::allocate()
       {
         auto size = static_cast<size_t>(size_);
-        f_.resize(size);
-        y_.resize(size);
-        yp_.resize(size);
+        this->allocateState(size);
         tag_.resize(size);
         variable_indices_.resize(size);
         residual_indices_.resize(size);
diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp
index 29c4b07bd..1f653849b 100644
--- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp
+++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp
@@ -136,9 +136,7 @@ namespace GridKit
       {
         // Allocate local component data
         auto size = static_cast<size_t>(size_); // avoid compiler warnings
-        f_.resize(size);
-        y_.resize(size);
-        yp_.resize(size);
+        this->allocateState(size);
         tag_.resize(size);
         variable_indices_.resize(size);
         residual_indices_.resize(size);
diff --git a/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp b/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp
index 58ea9c1f5..2d944b549 100644
--- a/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Load/LoadEnzyme.cpp
@@ -59,7 +59,7 @@ namespace GridKit
                                                  (bus_->getVariableIndices()).data(),
                                                  y_.data(),
                                                  yp_.data(),
-                                                 (bus_->y()).data(),
+                                                 bus_->yData(),
                                                  J_rows_buffer_,
                                                  J_cols_buffer_,
                                                  J_vals_buffer_,
diff --git a/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp
index e124bc375..5380969de 100644
--- a/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Load/LoadImpl.hpp
@@ -86,9 +86,7 @@ namespace GridKit
       // std::cout << "Allocate Load..." << std::endl;
 
       auto size = static_cast<size_t>(size_); // avoid compiler warnings
-      f_.resize(size);
-      y_.resize(size);
-      yp_.resize(size);
+      this->allocateState(size);
       tag_.resize(size);
       variable_indices_.resize(size);
       residual_indices_.resize(size);
diff --git a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp
index a5f7e6dfc..758753998 100644
--- a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPEnzyme.cpp
@@ -54,7 +54,7 @@ namespace GridKit
                                                  (bus_->getVariableIndices()).data(),
                                                  y_.data(),
                                                  yp_.data(),
-                                                 (bus_->y()).data(),
+                                                 bus_->yData(),
                                                  J_rows_buffer_,
                                                  J_cols_buffer_,
                                                  J_vals_buffer_,
diff --git a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPImpl.hpp b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPImpl.hpp
index 8dc38dc83..98911c996 100644
--- a/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPImpl.hpp
@@ -102,9 +102,7 @@ namespace GridKit
       // std::cout << "Allocate Load..." << std::endl;
 
       auto size = static_cast<size_t>(size_); // avoid compiler warnings
-      f_.resize(size);
-      y_.resize(size);
-      yp_.resize(size);
+      this->allocateState(size);
       tag_.resize(size);
       variable_indices_.resize(size);
       residual_indices_.resize(size);
diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp
index 70b49977d..ec7cc719e 100644
--- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp
@@ -157,9 +157,7 @@ namespace GridKit
       int Ieeest<scalar_type, index_type>::allocate()
       {
         auto size = static_cast<size_t>(size_);
-        f_.resize(size);
-        y_.resize(size);
-        yp_.resize(size);
+        this->allocateState(size);
         tag_.resize(size);
         variable_indices_.resize(size);
         residual_indices_.resize(size);
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp
index 2fb91a858..3190a22e0 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouEnzyme.cpp
@@ -95,7 +95,7 @@ namespace GridKit
                                                  (bus_->getVariableIndices()).data(),
                                                  y_.data(),
                                                  yp_.data(),
-                                                 (bus_->y()).data(),
+                                                 bus_->yData(),
                                                  J_rows_buffer_,
                                                  J_cols_buffer_,
                                                  J_vals_buffer_,
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp
index 9cbb0b587..125e12968 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp
@@ -327,9 +327,7 @@ namespace GridKit
     {
       // Resize component model data
       auto size = static_cast<size_t>(size_);
-      f_.resize(size);
-      y_.resize(size);
-      yp_.resize(size);
+      this->allocateState(size);
       tag_.resize(size);
       variable_indices_.resize(size);
       residual_indices_.resize(size);
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp
index deb4b1f5a..ae45672ad 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp
@@ -95,7 +95,7 @@ namespace GridKit
                                                  (bus_->getVariableIndices()).data(),
                                                  y_.data(),
                                                  yp_.data(),
-                                                 (bus_->y()).data(),
+                                                 bus_->yData(),
                                                  J_rows_buffer_,
                                                  J_cols_buffer_,
                                                  J_vals_buffer_,
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp
index 9b779ded5..a6243fb9d 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp
@@ -196,9 +196,7 @@ namespace GridKit
     {
       // Resize component model data
       auto size = static_cast<size_t>(size_);
-      f_.resize(size);
-      y_.resize(size);
-      yp_.resize(size);
+      this->allocateState(size);
       tag_.resize(size);
       variable_indices_.resize(size);
       residual_indices_.resize(size);
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp
index f5f664035..135dcf920 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalEnzyme.cpp
@@ -60,7 +60,7 @@ namespace GridKit
                                                  (bus_->getVariableIndices()).data(),
                                                  y_.data(),
                                                  yp_.data(),
-                                                 (bus_->y()).data(),
+                                                 bus_->yData(),
                                                  J_rows_buffer_,
                                                  J_cols_buffer_,
                                                  J_vals_buffer_,
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp
index 37578193a..307882703 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp
@@ -180,9 +180,7 @@ namespace GridKit
     {
       // Resize component model data
       auto size = static_cast<size_t>(size_);
-      f_.resize(size);
-      y_.resize(size);
-      yp_.resize(size);
+      this->allocateState(size);
       tag_.resize(size);
       variable_indices_.resize(size);
       residual_indices_.resize(size);
diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp
index d7b4cc9d7..3db97eeaf 100644
--- a/GridKit/Model/PhasorDynamics/SystemModel.hpp
+++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp
@@ -1,6 +1,7 @@
 #pragma once
 
 #include <cassert>
+#include <cstddef>
 #include <iostream>
 #include <vector>
 
@@ -48,13 +49,14 @@ namespace GridKit
       using PhasorDynamics::Component<scalar_type, index_type>::residual_indices_;
 
     public:
-      using ScalarT    = scalar_type;
-      using IdxT       = index_type;
-      using RealT      = typename Model::Evaluator<ScalarT, IdxT>::RealT;
-      using CsrMatrixT = typename Model::Evaluator<ScalarT, IdxT>::CsrMatrixT;
-      using BusT       = PhasorDynamics::BusBase<ScalarT, IdxT>;
-      using SignalT    = PhasorDynamics::SignalNode<ScalarT, IdxT>;
-      using ComponentT = PhasorDynamics::Component<ScalarT, IdxT>;
+      using ScalarT     = scalar_type;
+      using IdxT        = index_type;
+      using RealT       = typename Model::Evaluator<ScalarT, IdxT>::RealT;
+      using CsrMatrixT  = typename Model::Evaluator<ScalarT, IdxT>::CsrMatrixT;
+      using BusT        = PhasorDynamics::BusBase<ScalarT, IdxT>;
+      using SignalT     = PhasorDynamics::SignalNode<ScalarT, IdxT>;
+      using ComponentT  = PhasorDynamics::Component<ScalarT, IdxT>;
+      using VectorViewT = typename ComponentT::VectorViewT;
 
       /**
        * @brief Constructor for the system model
@@ -443,9 +445,9 @@ namespace GridKit
       /**
        * @brief Allocate buses, components, and system objects.
        *
-       * This method first allocates bus objects, then component objects,
-       * and computes system size (number of unknowns). Once the size is
-       * computed, system global objects are allocated.
+       * This method computes system size from bus/component sizes, allocates
+       * global objects, binds each element to its global slice, then allocates
+       * bus and component objects.
        *
        * @post size_ >= 1
        *
@@ -453,35 +455,16 @@ namespace GridKit
       int allocate() override
       {
         size_ = 0;
-
-        // Allocate all buses
         for (const auto& bus : buses_)
         {
-          bus->allocate();
-          for (IdxT j = 0; j < bus->size(); ++j)
-          {
-            bus->setVariableIndex(j, size_ + j);
-            bus->setResidualIndex(j, size_ + j);
-          }
           size_ += bus->size();
         }
-
-        // Allocate all components
         for (const auto& component : components_)
         {
-          component->allocate();
-          for (IdxT j = 0; j < component->size(); ++j)
-          {
-            component->setVariableIndex(j, size_ + j);
-            component->setResidualIndex(j, size_ + j);
-          }
           size_ += component->size();
         }
 
-        // Allocate global vectors
-        y_.resize(size_);
-        yp_.resize(size_);
-        f_.resize(size_);
+        this->allocateState(static_cast<std::size_t>(size_));
         tag_.resize(size_);
         variable_indices_.resize(size_);
         residual_indices_.resize(size_);
@@ -493,6 +476,41 @@ namespace GridKit
           this->setResidualIndex(j, j);
         }
 
+        auto bind = [this](auto* element, IdxT off, IdxT n)
+        {
+          const auto o    = static_cast<std::size_t>(off);
+          const auto size = static_cast<std::size_t>(n);
+          assert(o + size <= y_.size());
+          element->bindStateViews(VectorViewT(size == 0 ? nullptr : y_.data() + o, size),
+                                  VectorViewT(size == 0 ? nullptr : yp_.data() + o, size),
+                                  VectorViewT(size == 0 ? nullptr : f_.data() + o, size));
+        };
+        IdxT off = 0;
+        for (const auto& bus : buses_)
+        {
+          const IdxT n = bus->size();
+          bind(bus, off, n);
+          bus->allocate();
+          for (IdxT j = 0; j < n; ++j)
+          {
+            bus->setVariableIndex(j, off + j);
+            bus->setResidualIndex(j, off + j);
+          }
+          off += n;
+        }
+        for (const auto& component : components_)
+        {
+          const IdxT n = component->size();
+          bind(component, off, n);
+          component->allocate();
+          for (IdxT j = 0; j < n; ++j)
+          {
+            component->setVariableIndex(j, off + j);
+            component->setResidualIndex(j, off + j);
+          }
+          off += n;
+        }
+
         // Verify component configuration
         int errorCount = this->verify();
         if (errorCount > 0)
@@ -579,42 +597,24 @@ namespace GridKit
        * Also, generators may write to control devices (e.g. governors,
        * exciters, etc.) during the initialization.
        *
-       * @todo Implement writting to system vectors in a thread-safe way.
+       * @todo Make initialization writes to system vectors thread-safe.
        *
-       * @note Currently assuming each component stores variables contiguously in memory and
-       * that these are simply concateneted in the global system.
+       * @note Each bus and component is bound to a contiguous slice of the global system state.
        */
       int initialize() override
       {
+        // Buses initialize first: components may read bus values during their own
+        // initialization. Writes land directly in the global vector via the bound views.
         for (const auto& bus : buses_)
         {
           bus->initialize();
         }
 
-        for (const auto& bus : buses_)
-        {
-          for (IdxT j = 0; j < bus->size(); ++j)
-          {
-            y_[bus->getVariableIndex(j)]  = bus->y()[j];
-            yp_[bus->getVariableIndex(j)] = bus->yp()[j];
-          }
-        }
-
-        // Initialize components
         for (const auto& component : components_)
         {
           component->initialize();
         }
 
-        for (const auto& component : components_)
-        {
-          for (IdxT j = 0; j < component->size(); ++j)
-          {
-            y_[component->getVariableIndex(j)]  = component->y()[j];
-            yp_[component->getVariableIndex(j)] = component->yp()[j];
-          }
-        }
-
         return 0;
       }
 
@@ -686,24 +686,19 @@ namespace GridKit
       /**
        * @brief Compute system residual vector
        *
-       * First, update bus and component variables from the system solution
-       * vector. Next, evaluate residuals in buses and components, and
-       * then copy values to the global residual vector.
+       * Buses and components are bound to the system state during allocation.
+       * Residual evaluation therefore reads the current system values and
+       * writes residuals directly into the global residual vector.
        *
        * @warning Residuals must be computed for buses, before component
        * residuals are computed. Buses own residuals for currents
        * Ir and Ii, but the contributions to these residuals come
        * from components. Buses assign their residual values, while components
-       * add to those values by in-place adition. This is why (for now) bus
+       * add to those values by in-place addition. This is why (for now) bus
        * residuals need to be computed first.
-       *
-       * @todo Here, components write to local values, which are then copied
-       * to global system vectors. Make components write to the system
-       * vectors directly.
        */
       int evaluateResidual() override
       {
-        updateVariables();
 
         for (const auto& bus : buses_)
         {
@@ -715,23 +710,6 @@ namespace GridKit
           component->evaluateResidual();
         }
 
-        // Update residual vector
-        for (const auto& bus : buses_)
-        {
-          for (IdxT j = 0; j < bus->size(); ++j)
-          {
-            f_[bus->getResidualIndex(j)] = bus->getResidual()[j];
-          }
-        }
-
-        for (const auto& component : components_)
-        {
-          for (IdxT j = 0; j < component->size(); ++j)
-          {
-            f_[component->getResidualIndex(j)] = component->getResidual()[j];
-          }
-        }
-
         return 0;
       }
 
@@ -769,7 +747,7 @@ namespace GridKit
           IdxT nnz_dup = 0;
           for (const auto& component : components_)
           {
-            auto component_jacobian = component->getJacobian();
+            auto& component_jacobian = component->getJacobian();
 
             std::tuple<std::vector<IdxT>&, std::vector<IdxT>&, std::vector<RealT>&> component_jacobian_entries = component_jacobian.getEntries(false);
             const auto [rows, columns, values]                                                                 = component_jacobian_entries;
@@ -780,7 +758,7 @@ namespace GridKit
           }
           for (const auto& bus : buses_)
           {
-            auto bus_jacobian = bus->getJacobian();
+            auto& bus_jacobian = bus->getJacobian();
 
             std::tuple<std::vector<IdxT>&, std::vector<IdxT>&, std::vector<RealT>&> bus_jacobian_entries = bus_jacobian.getEntries(false);
             const auto [rows, columns, values]                                                           = bus_jacobian_entries;
@@ -798,7 +776,7 @@ namespace GridKit
           IdxT counter = 0;
           for (const auto& component : components_)
           {
-            auto component_jacobian = component->getJacobian();
+            auto& component_jacobian = component->getJacobian();
 
             std::tuple<std::vector<IdxT>&, std::vector<IdxT>&, std::vector<RealT>&> component_jacobian_entries = component_jacobian.getEntries(false);
             const auto [rows, columns, values]                                                                 = component_jacobian_entries;
@@ -812,7 +790,7 @@ namespace GridKit
           }
           for (const auto& bus : buses_)
           {
-            auto bus_jacobian = bus->getJacobian();
+            auto& bus_jacobian = bus->getJacobian();
 
             std::tuple<std::vector<IdxT>&, std::vector<IdxT>&, std::vector<RealT>&> bus_jacobian_entries = bus_jacobian.getEntries(false);
             const auto [rows, columns, values]                                                           = bus_jacobian_entries;
@@ -867,7 +845,7 @@ namespace GridKit
           IdxT counter = 0;
           for (const auto& component : components_)
           {
-            auto component_jacobian = component->getJacobian();
+            auto& component_jacobian = component->getJacobian();
 
             std::tuple<std::vector<IdxT>&, std::vector<IdxT>&, std::vector<RealT>&> component_jacobian_entries = component_jacobian.getEntries(false);
             const auto [rows, columns, values]                                                                 = component_jacobian_entries;
@@ -879,7 +857,7 @@ namespace GridKit
           }
           for (const auto& bus : buses_)
           {
-            auto bus_jacobian = bus->getJacobian();
+            auto& bus_jacobian = bus->getJacobian();
 
             std::tuple<std::vector<IdxT>&, std::vector<IdxT>&, std::vector<RealT>&> bus_jacobian_entries = bus_jacobian.getEntries(false);
             const auto [rows, columns, values]                                                           = bus_jacobian_entries;
@@ -912,29 +890,6 @@ namespace GridKit
         monitor_.print();
       }
 
-      /**
-       * @brief Update variables in buses and components
-       */
-      void updateVariables()
-      {
-        for (const auto& bus : buses_)
-        {
-          for (IdxT j = 0; j < bus->size(); ++j)
-          {
-            bus->y()[j]  = y_[bus->getVariableIndex(j)];
-            bus->yp()[j] = yp_[bus->getVariableIndex(j)];
-          }
-        }
-        for (const auto& component : components_)
-        {
-          for (IdxT j = 0; j < component->size(); ++j)
-          {
-            component->y()[j]  = y_[component->getVariableIndex(j)];
-            component->yp()[j] = yp_[component->getVariableIndex(j)];
-          }
-        }
-      }
-
       /**
        * @brief Update time
        *
@@ -947,8 +902,6 @@ namespace GridKit
         {
           component->updateTime(t, a);
         }
-
-        updateVariables();
       }
 
       /**
diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp
index 94bd68412..881fe97ea 100644
--- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp
+++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp
@@ -1,8 +1,10 @@
 #pragma once
 
+#include <cstddef>
 #include <iomanip>
 #include <iostream>
 #include <sstream>
+#include <stdexcept>
 #include <string>
 
 #include <GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp>
@@ -11,6 +13,7 @@
 #include <GridKit/Model/PhasorDynamics/Branch/BranchData.hpp>
 #include <GridKit/Model/PhasorDynamics/Bus/Bus.hpp>
 #include <GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp>
+#include <GridKit/Model/PhasorDynamics/Load/Load.hpp>
 #include <GridKit/Model/PhasorDynamics/SystemModel.hpp>
 #include <GridKit/Model/PhasorDynamics/SystemModelData.hpp>
 #include <GridKit/Testing/TestHelpers.hpp>
@@ -153,6 +156,54 @@ namespace GridKit
         return success.report(__func__);
       }
 
+      TestOutcome elementStateAliasesSystemState()
+      {
+        TestStatus success = true;
+
+        PhasorDynamics::SystemModel<ScalarT, IdxT> system;
+
+        PhasorDynamics::Bus<ScalarT, IdxT> bus1(1.0, 2.0);
+        PhasorDynamics::Bus<ScalarT, IdxT> bus2(3.0, 4.0);
+        system.addBus(&bus1);
+        system.addBus(&bus2);
+
+        PhasorDynamics::Branch<ScalarT, IdxT> branch(&bus1, &bus2, 0.0, 1.0, 0.0, 0.0);
+        system.addComponent(&branch);
+
+        PhasorDynamics::Load<ScalarT, IdxT> load(&bus2, 1.0, 1.0);
+        system.addComponent(&load);
+
+        system.allocate();
+        system.initialize();
+
+        const auto bus2_vr_index = static_cast<std::size_t>(bus2.getVariableIndex(0));
+        const auto bus2_vi_index = static_cast<std::size_t>(bus2.getVariableIndex(1));
+        const auto load_ir_index = static_cast<std::size_t>(load.getVariableIndex(0));
+        const auto load_ii_index = static_cast<std::size_t>(load.getVariableIndex(1));
+
+        system.y()[bus2_vr_index]  = 5.0;
+        system.y()[bus2_vi_index]  = 6.0;
+        success                   *= isEqual(bus2.Vr(), static_cast<ScalarT>(5.0));
+        success                   *= isEqual(bus2.Vi(), static_cast<ScalarT>(6.0));
+
+        bus2.Vr()  = 7.0;
+        bus2.Vi()  = 8.0;
+        success   *= isEqual(system.y()[bus2_vr_index], static_cast<ScalarT>(7.0));
+        success   *= isEqual(system.y()[bus2_vi_index], static_cast<ScalarT>(8.0));
+
+        system.y()[load_ir_index]  = 9.0;
+        system.y()[load_ii_index]  = 10.0;
+        success                   *= isEqual(load.yData()[0], static_cast<ScalarT>(9.0));
+        success                   *= isEqual(load.yData()[1], static_cast<ScalarT>(10.0));
+
+        load.yData()[0]  = 11.0;
+        load.yData()[1]  = 12.0;
+        success         *= isEqual(system.y()[load_ir_index], static_cast<ScalarT>(11.0));
+        success         *= isEqual(system.y()[load_ii_index], static_cast<ScalarT>(12.0));
+
+        return success.report(__func__);
+      }
+
       /**
        * @brief Test for exception when signals are incorrectly configured
        */
diff --git a/tests/UnitTests/PhasorDynamics/runSystemTests.cpp b/tests/UnitTests/PhasorDynamics/runSystemTests.cpp
index a8b432be2..b337028aa 100644
--- a/tests/UnitTests/PhasorDynamics/runSystemTests.cpp
+++ b/tests/UnitTests/PhasorDynamics/runSystemTests.cpp
@@ -10,6 +10,7 @@ int main()
 
   result += test.constructor();
   result += test.composer();
+  result += test.elementStateAliasesSystemState();
 #ifdef GRIDKIT_ENABLE_ENZYME
   result += test.jacobian();
 #endif