diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt
index 669af876e8da765ec3ac0d483a38c76b4e7ad414..fda3d1b2c485880a882f7948a84f0c2c6bbd697b 100644
--- a/src/analysis/CMakeLists.txt
+++ b/src/analysis/CMakeLists.txt
@@ -2,5 +2,6 @@
 
 add_library(
   PugsAnalysis
+  Polynomial1DRootComputer.cpp
   SturmSequence.cpp
 )
diff --git a/src/analysis/Polynomial1DRootComputer.cpp b/src/analysis/Polynomial1DRootComputer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c86d2ad3117a734979d41b9161cef3d7e8a4177e
--- /dev/null
+++ b/src/analysis/Polynomial1DRootComputer.cpp
@@ -0,0 +1,64 @@
+#include <analysis/Polynomial1DRootComputer.hpp>
+
+double
+Polynomial1DRootComputer::_computeRootBetween(const PositionSignChanges& left_bound,
+                                              const PositionSignChanges& right_bound) const
+{
+  double left  = left_bound.position;
+  double right = right_bound.position;
+
+  const double interval_size = right - left;
+
+  const Polynomial1D& p0 = m_sturm_sequence.p(0);
+  double left_value      = p0(left);
+
+  int cpt = 0;
+  while (right - left > 1E-12 * interval_size) {
+    double middle       = 0.5 * (right + left);
+    double middle_value = p0(middle);
+    if (middle_value * left_value < 0) {
+      right = middle;
+    } else {
+      left_value = middle_value;
+      left       = middle;
+    }
+    cpt++;
+  }
+
+  return 0.5 * (left + right);
+}
+
+std::optional<double>
+Polynomial1DRootComputer::getFirstRoot() const
+{
+  PositionSignChanges l = {m_left, m_sturm_sequence.numberOfSignChanges(m_left)};
+  PositionSignChanges r = {m_right, m_sturm_sequence.numberOfSignChanges(m_right)};
+
+  if (l.nb_sign_changes - r.nb_sign_changes > 0) {
+    if (l.nb_sign_changes - r.nb_sign_changes > 1) {
+      PositionSignChanges a = l;
+      PositionSignChanges b = r;
+      do {
+        Assert(a.nb_sign_changes - b.nb_sign_changes > 0);
+
+        double position_k      = 0.5 * (a.position + b.position);
+        PositionSignChanges rk = {position_k, m_sturm_sequence.numberOfSignChanges(position_k)};
+
+        if (l.nb_sign_changes - rk.nb_sign_changes > 0) {
+          b = rk;
+        } else if (l.nb_sign_changes - rk.nb_sign_changes == 0) {
+          a = rk;
+        }
+      } while (a.nb_sign_changes - b.nb_sign_changes != 1);
+      return _computeRootBetween(a, b);
+    } else {
+      return _computeRootBetween(l, r);
+    }
+  } else {
+    return {};
+  }
+}
+
+Polynomial1DRootComputer::Polynomial1DRootComputer(const Polynomial1D& p, const double a, const double b)
+  : m_sturm_sequence{p}, m_left{a}, m_right{b}
+{}
diff --git a/src/analysis/Polynomial1DRootComputer.hpp b/src/analysis/Polynomial1DRootComputer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0664216e58223bb49aacdb7314074b0067991a53
--- /dev/null
+++ b/src/analysis/Polynomial1DRootComputer.hpp
@@ -0,0 +1,31 @@
+#ifndef POLYNOMIAL_1D_ROOT_COMPUTER_HPP
+#define POLYNOMIAL_1D_ROOT_COMPUTER_HPP
+
+#include <analysis/SturmSequence.hpp>
+
+#include <optional>
+
+class Polynomial1DRootComputer
+{
+ private:
+  struct PositionSignChanges
+  {
+    double position;
+    size_t nb_sign_changes;
+  };
+
+  const SturmSequence m_sturm_sequence;
+  const double m_left;
+  const double m_right;
+
+  double _computeRootBetween(const PositionSignChanges& l, const PositionSignChanges& r) const;
+
+ public:
+  std::optional<double> getFirstRoot() const;
+
+  Polynomial1DRootComputer(const Polynomial1D& p, const double a, const double b);
+
+  ~Polynomial1DRootComputer() = default;
+};
+
+#endif   // POLYNOMIAL_1D_ROOT_COMPUTER_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 68f1bba36dcc0aa51427c719be2259e959788dfa..4a602a294bf663d41e075faf29ebaa943058ca54 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -84,6 +84,7 @@ add_executable (unit_tests
   test_OStreamProcessor.cpp
   test_ParseError.cpp
   test_Polynomial1D.cpp
+  test_Polynomial1DRootComputer.cpp
   test_PugsAssert.cpp
   test_PugsFunctionAdapter.cpp
   test_PugsUtils.cpp
diff --git a/tests/test_Polynomial1DRootComputer.cpp b/tests/test_Polynomial1DRootComputer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e17509f29596c63705ddaa0b313b21d5aeffb24c
--- /dev/null
+++ b/tests/test_Polynomial1DRootComputer.cpp
@@ -0,0 +1,45 @@
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include <analysis/Polynomial1DRootComputer.hpp>
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("Polynomial1DRootComputer", "[analysis]")
+{
+  SECTION("No root in interval")
+  {
+    Polynomial1D p({-1, 0, 0, 1});
+
+    Polynomial1DRootComputer computer(p, -4, -3);
+
+    std::optional optional_root = computer.getFirstRoot();
+
+    REQUIRE(optional_root.has_value() == false);
+  }
+
+  SECTION("one root in interval")
+  {
+    Polynomial1D p({-1, 0, 0, 1});
+
+    Polynomial1DRootComputer computer(p, -0.3, 1.2);
+
+    std::optional optional_root = computer.getFirstRoot();
+
+    REQUIRE(optional_root.has_value() == true);
+    REQUIRE(optional_root.value() == Catch::Approx(1));
+  }
+
+  SECTION("multiple roots in interval")
+  {
+    Polynomial1D p =
+      Polynomial1D{{-0.5, 1}} * Polynomial1D{{-0.23, 1}} * Polynomial1D{{-0.27, 1}} * Polynomial1D{{-0.43, 1}};
+
+    Polynomial1DRootComputer computer(p, -0.3, 1.2);
+
+    std::optional optional_root = computer.getFirstRoot();
+
+    REQUIRE(optional_root.has_value() == true);
+    REQUIRE(optional_root.value() == Catch::Approx(0.23));
+  }
+}