diff --git a/CMakeLists.txt b/CMakeLists.txt
index a4aa863b266890293cb48ceaf21990120f4d303b..67136d4686b0caa3063ec9090af053ebd46ee136 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -545,6 +545,7 @@ target_link_libraries(
   pugs
   PugsMesh
   PugsAlgebra
+  PugsAnalysis
   PugsUtils
   PugsLanguage
   PugsLanguageAST
@@ -552,6 +553,7 @@ target_link_libraries(
   PugsLanguageAlgorithms
   PugsMesh
   PugsAlgebra
+  PugsAnalysis
   PugsScheme
   PugsUtils
   PugsOutput
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 16904a7be0cddc9aa54206638ab4b5c759cd0de8..54cae625bc7c899c433bfdf4021ef6242dd5aaf4 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -12,6 +12,9 @@ add_subdirectory(language)
 # Pugs algebra
 add_subdirectory(algebra)
 
+# Pugs analysis
+add_subdirectory(analysis)
+
 # Pugs mesh
 add_subdirectory(mesh)
 
diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..669af876e8da765ec3ac0d483a38c76b4e7ad414
--- /dev/null
+++ b/src/analysis/CMakeLists.txt
@@ -0,0 +1,6 @@
+# ------------------- Source files --------------------
+
+add_library(
+  PugsAnalysis
+  SturmSequence.cpp
+)
diff --git a/src/analysis/SturmSequence.cpp b/src/analysis/SturmSequence.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bda7f5fa58be9bfcee2c2ef1b207583053ba8eb5
--- /dev/null
+++ b/src/analysis/SturmSequence.cpp
@@ -0,0 +1,39 @@
+#include <analysis/SturmSequence.hpp>
+
+SturmSequence::SturmSequence(const Polynomial1D& p0) : SturmSequence{Polynomial1D{p0}} {}
+
+SturmSequence::SturmSequence(Polynomial1D&& p0)
+{
+  m_sequence.reserve(p0.degree());
+  m_sequence.emplace_back(p0);
+  Polynomial1D p1 = derive(m_sequence[0]);
+  if (p1.degree() > 0 or (std::abs(p1.coefficient(0)) > 1E-14)) {
+    m_sequence.emplace_back(std::move(p1));
+
+    do {
+      Polynomial1D pk = -(m_sequence[m_sequence.size() - 2] % m_sequence[m_sequence.size() - 1]);
+      if (pk.degree() > 0 or (std::abs(pk.coefficient(0)) > 1E-14)) {
+        m_sequence.emplace_back(std::move(pk));
+      } else {
+        break;
+      }
+    } while (true);
+  }
+}
+
+size_t
+SturmSequence::numberOfSignChanges(const double x)
+{
+  std::vector<double> values;
+  values.reserve(m_sequence.size());
+  for (const auto& pk : m_sequence) {
+    values.emplace_back(pk(x));
+  }
+
+  size_t number_of_sign_changes = 0;
+  for (size_t i = 1; i < m_sequence.size(); ++i) {
+    number_of_sign_changes += (values[i - 1] * values[i] < 0);
+  }
+
+  return number_of_sign_changes;
+}
diff --git a/src/analysis/SturmSequence.hpp b/src/analysis/SturmSequence.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5445cb65f96887faca185f713e8cac43c93fc312
--- /dev/null
+++ b/src/analysis/SturmSequence.hpp
@@ -0,0 +1,32 @@
+#ifndef STURM_SEQUENCE_HPP
+#define STURM_SEQUENCE_HPP
+
+#include <analysis/Polynomial1D.hpp>
+
+#include <vector>
+
+class SturmSequence
+{
+ private:
+  std::vector<Polynomial1D> m_sequence;
+
+ public:
+  PUGS_INLINE
+  size_t
+  size() const
+  {
+    return m_sequence.size();
+  }
+
+  size_t numberOfSignChanges(const double x);
+
+  SturmSequence(const Polynomial1D& p0);
+  SturmSequence(Polynomial1D&& p0);
+
+  SturmSequence(SturmSequence&&)      = delete;
+  SturmSequence(const SturmSequence&) = delete;
+
+  ~SturmSequence() = default;
+};
+
+#endif   // STURM_SEQUENCE_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 827da322463b0b94ef9ab0ef1aa6714211c53c27..68f1bba36dcc0aa51427c719be2259e959788dfa 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -89,6 +89,7 @@ add_executable (unit_tests
   test_PugsUtils.cpp
   test_RevisionInfo.cpp
   test_SparseMatrixDescriptor.cpp
+  test_SturmSequence.cpp
   test_SymbolTable.cpp
   test_Table.cpp
   test_Timer.cpp
@@ -128,6 +129,7 @@ target_link_libraries (unit_tests
   PugsLanguage
   PugsMesh
   PugsAlgebra
+  PugsAnalysis
   PugsScheme
   PugsOutput
   PugsUtils
@@ -150,6 +152,7 @@ target_link_libraries (mpi_unit_tests
   PugsLanguageAlgorithms
   PugsMesh
   PugsAlgebra
+  PugsAnalysis
   PugsUtils
   PugsLanguageUtils
   PugsScheme
diff --git a/tests/test_SturmSequence.cpp b/tests/test_SturmSequence.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..198fd029f0cb0f80e9f0dff03aaf7068723f3065
--- /dev/null
+++ b/tests/test_SturmSequence.cpp
@@ -0,0 +1,18 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/SturmSequence.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("SturmSequence", "[analysis]")
+{
+  Polynomial1D p({-1, -1, 0, 1, 1});
+
+  SturmSequence sturm_sequence(p);
+
+  REQUIRE(sturm_sequence.size() == 5);
+
+  REQUIRE(sturm_sequence.numberOfSignChanges(-10) == 3);
+  REQUIRE(sturm_sequence.numberOfSignChanges(10) == 1);
+}