diff --git a/src/LinearAlgebra/SparsityPattern/Variable.hpp b/src/LinearAlgebra/SparsityPattern/Variable.hpp
index 4015f65a..5b95eb61 100644
--- a/src/LinearAlgebra/SparsityPattern/Variable.hpp
+++ b/src/LinearAlgebra/SparsityPattern/Variable.hpp
@@ -14,6 +14,8 @@
 #include <string>
 #include <vector>
 
+#include <ScalarTraits.hpp>
+
 namespace GridKit
 {
   namespace Sparse
@@ -326,5 +328,18 @@ namespace GridKit
   } // namespace Sparse
 } // namespace GridKit
 
+namespace GridKit
+{
+  template <>
+  class ScalarTraits<Sparse::Variable>
+  {
+  public:
+    typedef double real_type;
+    typedef double norm_type;
+    typedef double scalar_type;
+  };
+
+} // namespace GridKit
+
 #include "VariableImplementation.hpp"
 #include "VariableOperators.hpp"
diff --git a/src/Model/PhasorDynamics/Branch/Branch.cpp b/src/Model/PhasorDynamics/Branch/Branch.cpp
index 75ccb8d7..0ccbe1a1 100644
--- a/src/Model/PhasorDynamics/Branch/Branch.cpp
+++ b/src/Model/PhasorDynamics/Branch/Branch.cpp
@@ -221,6 +221,8 @@ namespace GridKit
     // Available template instantiations
     template class Branch<double, long int>;
     template class Branch<double, size_t>;
+    template class Branch<Sparse::Variable, long int>;
+    template class Branch<Sparse::Variable, size_t>;
 
   } // namespace PhasorDynamics
 } // namespace GridKit
diff --git a/src/Model/PhasorDynamics/Bus/Bus.cpp b/src/Model/PhasorDynamics/Bus/Bus.cpp
index f8005b2b..7253ee66 100644
--- a/src/Model/PhasorDynamics/Bus/Bus.cpp
+++ b/src/Model/PhasorDynamics/Bus/Bus.cpp
@@ -215,6 +215,8 @@ namespace GridKit
     // Available template instantiations
     template class Bus<double, long int>;
     template class Bus<double, size_t>;
+    template class Bus<Sparse::Variable, long int>;
+    template class Bus<Sparse::Variable, size_t>;
 
   } // namespace PhasorDynamics
 } // namespace GridKit
diff --git a/src/Model/PhasorDynamics/Bus/BusInfinite.cpp b/src/Model/PhasorDynamics/Bus/BusInfinite.cpp
index 6206813d..f4f18b83 100644
--- a/src/Model/PhasorDynamics/Bus/BusInfinite.cpp
+++ b/src/Model/PhasorDynamics/Bus/BusInfinite.cpp
@@ -195,6 +195,8 @@ namespace GridKit
     // Available template instantiations
     template class BusInfinite<double, long int>;
     template class BusInfinite<double, size_t>;
+    template class BusInfinite<Sparse::Variable, long int>;
+    template class BusInfinite<Sparse::Variable, size_t>;
 
   } // namespace PhasorDynamics
 } // namespace GridKit
diff --git a/src/Model/PhasorDynamics/BusBase.hpp b/src/Model/PhasorDynamics/BusBase.hpp
index 804cbfca..cd92d1bc 100644
--- a/src/Model/PhasorDynamics/BusBase.hpp
+++ b/src/Model/PhasorDynamics/BusBase.hpp
@@ -2,6 +2,7 @@
 
 #include <vector>
 
+#include <LinearAlgebra/SparsityPattern/Variable.hpp>
 #include <Model/Evaluator.hpp>
 
 namespace GridKit
diff --git a/src/Model/PhasorDynamics/Component.hpp b/src/Model/PhasorDynamics/Component.hpp
index c061ff62..eec865bf 100644
--- a/src/Model/PhasorDynamics/Component.hpp
+++ b/src/Model/PhasorDynamics/Component.hpp
@@ -2,6 +2,7 @@
 
 #include <vector>
 
+#include <LinearAlgebra/SparsityPattern/Variable.hpp>
 #include <Model/Evaluator.hpp>
 
 namespace GridKit
diff --git a/src/Model/PhasorDynamics/Load/Load.cpp b/src/Model/PhasorDynamics/Load/Load.cpp
index 84763e3d..8c042a84 100644
--- a/src/Model/PhasorDynamics/Load/Load.cpp
+++ b/src/Model/PhasorDynamics/Load/Load.cpp
@@ -173,6 +173,8 @@ namespace GridKit
     // Available template instantiations
     template class Load<double, long int>;
     template class Load<double, size_t>;
+    template class Load<Sparse::Variable, long int>;
+    template class Load<Sparse::Variable, size_t>;
 
   } // namespace PhasorDynamics
 } // namespace GridKit
diff --git a/src/Utilities/Testing.hpp b/src/Utilities/Testing.hpp
index ce57d2c5..6f91ef1f 100644
--- a/src/Utilities/Testing.hpp
+++ b/src/Utilities/Testing.hpp
@@ -10,6 +10,8 @@
 #include <cmath>
 #include <iostream>
 #include <limits>
+#include <map>
+#include <vector>
 
 #include <PowerSystemData.hpp>
 
@@ -220,6 +222,57 @@ namespace GridKit
       return fail == 0;
     }
 
+    /**
+     * @brief Equatlity comparison between maps with a tolerance for the scalar value
+     *
+     * @tparam IdxT
+     * @tparam RealT
+     *
+     * @param[in] a - first map to compare
+     * @param[in] b - second map to compare
+     * @param[in] tol - tolerance
+     * @return bool - true if the maps are equal; false otherwise
+     */
+    template <typename IdxT = size_t, typename RealT = double>
+    inline bool isEqual(std::map<IdxT, RealT> a,
+                        std::map<IdxT, RealT> b,
+                        const RealT           tol = std::numeric_limits<RealT>::epsilon())
+    {
+      int fail = 0;
+
+      if (a.size() != b.size())
+      {
+        fail++;
+        errs() << "Containers do not have the same size! "
+               << "a.size() = " << a.size() << ", and "
+               << "b.size() = " << b.size() << "\n";
+      }
+
+      for (const auto& pair_a : a)
+      {
+        auto it = b.find(pair_a.first);
+        if (it != b.end())
+        {
+          if (!isEqual(pair_a.second, it->second, tol))
+          {
+            fail++;
+            errs() << "Mismatching map values! "
+                   << "a.first = " << pair_a.first << ", "
+                   << "a.second = " << pair_a.second << ", and "
+                   << "b.second = " << it->second << "\n";
+          }
+        }
+        else
+        {
+          fail++;
+          errs() << "Entry not found in the second container! "
+                 << "a.first = " << pair_a.first << "\n";
+        }
+      }
+
+      return fail == 0;
+    }
+
   } // namespace Testing
 
 } // namespace GridKit
diff --git a/tests/UnitTests/PhasorDynamics/BranchTests.hpp b/tests/UnitTests/PhasorDynamics/BranchTests.hpp
index 81e2c48e..2ca571e7 100644
--- a/tests/UnitTests/PhasorDynamics/BranchTests.hpp
+++ b/tests/UnitTests/PhasorDynamics/BranchTests.hpp
@@ -1,6 +1,7 @@
 #include <iomanip>
 #include <iostream>
 
+#include <LinearAlgebra/SparsityPattern/Variable.hpp>
 #include <Model/PhasorDynamics/Branch/Branch.hpp>
 #include <Model/PhasorDynamics/Bus/Bus.hpp>
 #include <Model/PhasorDynamics/Bus/BusInfinite.hpp>
@@ -77,6 +78,46 @@ namespace GridKit
         return success.report(__func__);
       }
 
+      TestOutcome jacobian()
+      {
+        TestStatus success = true;
+
+        real_type R{2.0}; ///< Branch series resistance
+        real_type X{4.0}; ///< Branch series reactance
+        real_type G{0.2}; ///< Branch shunt conductance
+        real_type B{1.2}; ///< Branch shunt charging
+
+        Sparse::Variable Vr1{10.0}; ///< Bus-1 real voltage
+        Sparse::Variable Vi1{20.0}; ///< Bus-1 imaginary voltage
+        Sparse::Variable Vr2{30.0}; ///< Bus-2 real voltage
+        Sparse::Variable Vi2{40.0}; ///< Bus-2 imaginary voltage
+
+        Vr1.setVariableNumber(0); ///< Independent variables: first
+        Vi1.setVariableNumber(1); ///< Independent variables: second
+        Vr2.setVariableNumber(2); ///< Independent variables: third
+        Vi2.setVariableNumber(3); ///< Independent variables: fourth
+
+        PhasorDynamics::BusInfinite<Sparse::Variable, IdxT> bus1(Vr1, Vi1);
+        PhasorDynamics::BusInfinite<Sparse::Variable, IdxT> bus2(Vr2, Vi2);
+
+        PhasorDynamics::Branch<Sparse::Variable, IdxT> branch(&bus1, &bus2, R, X, G, B);
+        branch.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking
+                                   ///< the dependencies
+
+        std::vector<Sparse::Variable>                residuals{bus1.Ir(), bus1.Ii(), bus2.Ir(), bus2.Ii()};
+        std::vector<Sparse::Variable::DependencyMap> ref = analyticalJacobian(R, X, G, B);
+
+        /// Compare dependencies computed automatically to the ones computed analytically
+        for (size_t i = 0; i < residuals.size(); ++i)
+        {
+          Sparse::Variable                       res           = residuals[i];
+          const Sparse::Variable::DependencyMap& dependencies  = res.getDependencies();
+          success                                             *= (GridKit::Testing::isEqual(dependencies, ref[i]));
+        }
+
+        return success.report(__func__);
+      }
+
       TestOutcome accessors()
       {
         TestStatus success = true;
@@ -139,6 +180,44 @@ namespace GridKit
 
         return success.report(__func__);
       }
+
+    private:
+      std::vector<Sparse::Variable::DependencyMap> analyticalJacobian(const real_type R,
+                                                                      const real_type X,
+                                                                      const real_type G,
+                                                                      const real_type B)
+      {
+        const real_type b = -X / (R * R + X * X);
+        const real_type g = R / (R * R + X * X);
+
+        real_type dIr1_dVr1 = -(g + 0.5 * G);
+        real_type dIr1_dVi1 = (b + 0.5 * B);
+        real_type dIr1_dVr2 = g;
+        real_type dIr1_dVi2 = -b;
+
+        real_type dIi1_dVr1 = -(b + 0.5 * B);
+        real_type dIi1_dVi1 = -(g + 0.5 * G);
+        real_type dIi1_dVr2 = b;
+        real_type dIi1_dVi2 = g;
+
+        real_type dIr2_dVr1 = g;
+        real_type dIr2_dVi1 = -b;
+        real_type dIr2_dVr2 = -(g + 0.5 * G);
+        real_type dIr2_dVi2 = (b + 0.5 * B);
+
+        real_type dIi2_dVr1 = b;
+        real_type dIi2_dVi1 = g;
+        real_type dIi2_dVr2 = -(b + 0.5 * B);
+        real_type dIi2_dVi2 = -(g + 0.5 * G);
+
+        std::vector<Sparse::Variable::DependencyMap> dependencies(4);
+        dependencies[0] = {{0, dIr1_dVr1}, {1, dIr1_dVi1}, {2, dIr1_dVr2}, {3, dIr1_dVi2}};
+        dependencies[1] = {{0, dIi1_dVr1}, {1, dIi1_dVi1}, {2, dIi1_dVr2}, {3, dIi1_dVi2}};
+        dependencies[2] = {{0, dIr2_dVr1}, {1, dIr2_dVi1}, {2, dIr2_dVr2}, {3, dIr2_dVi2}};
+        dependencies[3] = {{0, dIi2_dVr1}, {1, dIi2_dVi1}, {2, dIi2_dVr2}, {3, dIi2_dVi2}};
+
+        return dependencies;
+      }
     }; // class BranchTest
 
   } // namespace Testing
diff --git a/tests/UnitTests/PhasorDynamics/LoadTests.hpp b/tests/UnitTests/PhasorDynamics/LoadTests.hpp
index 7a9688bc..22c46aaf 100644
--- a/tests/UnitTests/PhasorDynamics/LoadTests.hpp
+++ b/tests/UnitTests/PhasorDynamics/LoadTests.hpp
@@ -3,6 +3,7 @@
 #include <iomanip>
 #include <iostream>
 
+#include <LinearAlgebra/SparsityPattern/Variable.hpp>
 #include <Model/PhasorDynamics/Bus/Bus.hpp>
 #include <Model/PhasorDynamics/Bus/BusInfinite.hpp>
 #include <Model/PhasorDynamics/Load/Load.hpp>
@@ -65,6 +66,59 @@ namespace GridKit
 
         return success.report(__func__);
       }
+
+      TestOutcome jacobian()
+      {
+        TestStatus success = true;
+
+        real_type R{2.0}; ///< Load resistance
+        real_type X{4.0}; ///< Load reactance
+
+        Sparse::Variable Vr{10.0}; ///< Bus real voltage
+        Sparse::Variable Vi{20.0}; ///< Bus imaginary voltage
+
+        Vr.setVariableNumber(0); ///< Independent variables: first
+        Vi.setVariableNumber(1); ///< Independent variables: second
+
+        PhasorDynamics::BusInfinite<Sparse::Variable, IdxT> bus(Vr, Vi);
+
+        PhasorDynamics::Load<Sparse::Variable, IdxT> load(&bus, R, X);
+        load.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking
+                                 ///< the dependencies
+
+        std::vector<Sparse::Variable>                residuals{bus.Ir(), bus.Ii()};
+        std::vector<Sparse::Variable::DependencyMap> ref = analyticalJacobian(R, X);
+
+        /// Compare dependencies computed automatically to the ones computed analytically
+        for (size_t i = 0; i < residuals.size(); ++i)
+        {
+          Sparse::Variable                       res           = residuals[i];
+          const Sparse::Variable::DependencyMap& dependencies  = res.getDependencies();
+          success                                             *= (GridKit::Testing::isEqual(dependencies, ref[i]));
+        }
+
+        return success.report(__func__);
+      }
+
+    private:
+      std::vector<Sparse::Variable::DependencyMap> analyticalJacobian(const real_type R,
+                                                                      const real_type X)
+      {
+        const real_type b = -X / (R * R + X * X);
+        const real_type g = R / (R * R + X * X);
+
+        real_type dIr_dVr = -g;
+        real_type dIr_dVi = -b;
+
+        real_type dIi_dVr = b;
+        real_type dIi_dVi = -g;
+
+        std::vector<Sparse::Variable::DependencyMap> dependencies(2);
+        dependencies[0] = {{0, dIr_dVr}, {1, dIr_dVi}};
+        dependencies[1] = {{0, dIi_dVr}, {1, dIi_dVi}};
+
+        return dependencies;
+      }
     };
 
   } // namespace Testing
diff --git a/tests/UnitTests/PhasorDynamics/runBranchTests.cpp b/tests/UnitTests/PhasorDynamics/runBranchTests.cpp
index b70a39ff..123f5d29 100644
--- a/tests/UnitTests/PhasorDynamics/runBranchTests.cpp
+++ b/tests/UnitTests/PhasorDynamics/runBranchTests.cpp
@@ -11,6 +11,7 @@ int main()
   result += test.constructor();
   result += test.accessors();
   result += test.residual();
+  result += test.jacobian();
 
   return result.summary();
 }
diff --git a/tests/UnitTests/PhasorDynamics/runLoadTests.cpp b/tests/UnitTests/PhasorDynamics/runLoadTests.cpp
index 7ad3daf0..41426a82 100644
--- a/tests/UnitTests/PhasorDynamics/runLoadTests.cpp
+++ b/tests/UnitTests/PhasorDynamics/runLoadTests.cpp
@@ -10,6 +10,7 @@ int main()
 
   result += test.constructor();
   result += test.residual();
+  result += test.jacobian();
 
   return result.summary();
 }