diff --git a/doc/interpolators/quintic_hermite.qbk b/doc/interpolators/quintic_hermite.qbk
index dafea2810..a2d550d4e 100644
--- a/doc/interpolators/quintic_hermite.qbk
+++ b/doc/interpolators/quintic_hermite.qbk
@@ -25,6 +25,8 @@ public:
 
     inline Real double_prime(Real x) const;
 
+    std::pair<Real, Real> domain() const;
+
     friend std::ostream& operator<<(std::ostream & os, const quintic_hermite & m);
 
     void push_back(Real x, Real y, Real dydx, Real d2ydx2);
@@ -41,6 +43,8 @@ public:
     inline Real prime(Real x) const;
 
     inline Real double_prime(Real x) const;
+
+    std::pair<Real, Real> domain() const;
 };
 
 template<class RandomAccessContainer>
@@ -56,6 +60,8 @@ public:
 
     inline Real double_prime(Real x) const;
 
+    std::pair<Real, Real> domain() const;
+
 }
 ``
 
diff --git a/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp b/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp
index 26a4e7498..28f04daca 100644
--- a/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp
+++ b/include/boost/math/interpolators/detail/quintic_hermite_detail.hpp
@@ -188,6 +188,11 @@ public:
         return 4*x_.size()*sizeof(x_);
     }
 
+    std::pair<Real, Real> domain() const
+    {
+        return {x_.front(), x_.back()};
+    }
+
 private:
     RandomAccessContainer x_;
     RandomAccessContainer y_;
@@ -257,7 +262,10 @@ public:
         Real ii = floor(s);
         auto i = static_cast<decltype(y_.size())>(ii);
         Real t = s - ii;
-
+        if (t == 0)
+        {
+            return y_[i];
+        }
         Real y0 = y_[i];
         Real y1 = y_[i+1];
         Real dy0 = dy_[i];
@@ -301,7 +309,10 @@ public:
         Real ii = floor(s);
         auto i = static_cast<decltype(y_.size())>(ii);
         Real t = s - ii;
-
+        if (t == 0)
+        {
+            return dy_[i]*inv_dx_;
+        }
         Real y0 = y_[i];
         Real y1 = y_[i+1];
         Real dy0 = dy_[i];
@@ -341,7 +352,10 @@ public:
         Real ii = floor(s);
         auto i = static_cast<decltype(y_.size())>(ii);
         Real t = s - ii;
-
+        if (t==0)
+        {
+            return d2y_[i]*2*inv_dx_*inv_dx_;
+        }
 
         Real y0 = y_[i];
         Real y1 = y_[i+1];
@@ -361,6 +375,12 @@ public:
         return 3*y_.size()*sizeof(Real) + 2*sizeof(Real);
     }
 
+    std::pair<Real, Real> domain() const
+    {
+        Real xf = x0_ + (y_.size()-1)/inv_dx_;
+        return {x0_, xf};
+    }
+
 private:
     RandomAccessContainer y_;
     RandomAccessContainer dy_;
@@ -425,6 +445,10 @@ public:
         Real ii = floor(s);
         auto i = static_cast<decltype(data_.size())>(ii);
         Real t = s - ii;
+        if (t == 0)
+        {
+            return data_[i][0];
+        }
 
         Real y0 = data_[i][0];
         Real dy0 = data_[i][1];
@@ -466,6 +490,10 @@ public:
         Real ii = floor(s);
         auto i = static_cast<decltype(data_.size())>(ii);
         Real t = s - ii;
+        if (t == 0)
+        {
+            return data_[i][1]*inv_dx_;
+        }
 
 
         Real y0 = data_[i][0];
@@ -507,8 +535,9 @@ public:
         Real ii = floor(s);
         auto i = static_cast<decltype(data_.size())>(ii);
         Real t = s - ii;
-
-
+        if (t == 0) {
+            return data_[i][2]*2*inv_dx_*inv_dx_;
+        }
         Real y0 = data_[i][0];
         Real dy0 = data_[i][1];
         Real d2y0 = data_[i][2];
@@ -527,6 +556,12 @@ public:
         return data_.size()*data_[0].size()*sizeof(Real) + 2*sizeof(Real);
     }
 
+    std::pair<Real, Real> domain() const
+    {
+        Real xf = x0_ + (data_.size()-1)/inv_dx_;
+        return {x0_, xf};
+    }
+
 private:
     RandomAccessContainer data_;
     Real x0_;
diff --git a/include/boost/math/interpolators/quintic_hermite.hpp b/include/boost/math/interpolators/quintic_hermite.hpp
index beb7238f1..627f02d28 100644
--- a/include/boost/math/interpolators/quintic_hermite.hpp
+++ b/include/boost/math/interpolators/quintic_hermite.hpp
@@ -46,6 +46,11 @@ public:
         return impl_->bytes() + sizeof(impl_);
     }
 
+    std::pair<Real, Real> domain() const
+    {
+        return impl_->domain();
+    }
+
 private:
     std::shared_ptr<detail::quintic_hermite_detail<RandomAccessContainer>> impl_;
 };
@@ -76,6 +81,11 @@ public:
         return impl_->bytes() + sizeof(impl_);
     }
 
+    std::pair<Real, Real> domain() const
+    {
+        return impl_->domain();
+    }
+
 private:
     std::shared_ptr<detail::cardinal_quintic_hermite_detail<RandomAccessContainer>> impl_;
 };
@@ -108,6 +118,11 @@ public:
     {
         return impl_->bytes() + sizeof(impl_);
     }
+
+    std::pair<Real, Real> domain() const
+    {
+        return impl_->domain();
+    }
 private:
     std::shared_ptr<detail::cardinal_quintic_hermite_detail_aos<RandomAccessContainer>> impl_;
 };
diff --git a/test/cubic_hermite_test.cpp b/test/cubic_hermite_test.cpp
index 4c9cfc003..0a23fb166 100644
--- a/test/cubic_hermite_test.cpp
+++ b/test/cubic_hermite_test.cpp
@@ -246,7 +246,7 @@ void test_cardinal_constant()
         }
     }
 
-        // Now check the boundaries:
+    // Now check the boundaries:
     Real tlo = x0;
     Real thi = x0 + (25-1)*dx;
     int samples = 5000;
diff --git a/test/quintic_hermite_test.cpp b/test/quintic_hermite_test.cpp
index a668e42e1..05299bae3 100644
--- a/test/quintic_hermite_test.cpp
+++ b/test/quintic_hermite_test.cpp
@@ -13,6 +13,7 @@
 #include <boost/random/uniform_real.hpp>
 #include <boost/random/mersenne_twister.hpp>
 #include <boost/math/interpolators/quintic_hermite.hpp>
+#include <boost/math/special_functions/next.hpp>
 #include <boost/circular_buffer.hpp>
 #ifdef BOOST_HAS_FLOAT128
 #include <boost/multiprecision/float128.hpp>
@@ -287,6 +288,25 @@ void test_cardinal_constant()
         CHECK_ULP_CLOSE(Real(0), qh_aos.prime(t), 24);
         CHECK_ULP_CLOSE(Real(0), qh_aos.double_prime(t), 24);
     }
+
+    // Now check the boundaries:
+    auto [tlo, thi] = qh.domain();
+    int samples = 5000;
+    int i = 0;
+    while (i++ < samples)
+    {
+        CHECK_ULP_CLOSE(Real(7), qh(tlo), 2);
+        CHECK_ULP_CLOSE(Real(7), qh(thi), 2);
+        CHECK_ULP_CLOSE(Real(7), qh_aos(tlo), 2);
+        CHECK_ULP_CLOSE(Real(7), qh_aos(thi), 2);
+        CHECK_ULP_CLOSE(Real(0), qh.prime(tlo), 2);
+        CHECK_ULP_CLOSE(Real(0), qh.prime(thi), 2);
+        CHECK_ULP_CLOSE(Real(0), qh_aos.prime(tlo), 2);
+        CHECK_ULP_CLOSE(Real(0), qh_aos.prime(thi), 2);
+
+        tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
+        thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
+    }
 }
 
 
@@ -322,6 +342,24 @@ void test_cardinal_linear()
         CHECK_ULP_CLOSE(Real(0), qh_aos.double_prime(t), 2);
     }
 
+    // Now check the boundaries:
+    auto [tlo, thi] = qh.domain();
+    int samples = 5000;
+    int i = 0;
+    while (i++ < samples)
+    {
+        CHECK_ULP_CLOSE(Real(tlo), qh(tlo), 2);
+        CHECK_ULP_CLOSE(Real(thi), qh(thi), 2);
+        CHECK_ULP_CLOSE(Real(tlo), qh_aos(tlo), 2);
+        CHECK_ULP_CLOSE(Real(thi), qh_aos(thi), 2);
+        CHECK_ULP_CLOSE(Real(1), qh.prime(tlo), 2);
+        CHECK_ULP_CLOSE(Real(1), qh.prime(thi), 128);
+        CHECK_ULP_CLOSE(Real(1), qh_aos.prime(tlo), 2);
+        CHECK_ULP_CLOSE(Real(1), qh_aos.prime(thi), 128);
+
+        tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
+        thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
+    }
 }
 
 template<typename Real>
@@ -366,6 +404,25 @@ void test_cardinal_quadratic()
         CHECK_ULP_CLOSE(t, qh_aos.prime(t), 12);
         CHECK_ULP_CLOSE(Real(1), qh_aos.double_prime(t), 64);
     }
+
+        // Now check the boundaries:
+    auto [tlo, thi] = qh.domain();
+    int samples = 5000;
+    int i = 0;
+    while (i++ < samples)
+    {
+        CHECK_ULP_CLOSE(tlo*tlo/2, qh(tlo), 16);
+        CHECK_ULP_CLOSE(thi*thi/2, qh(thi), 16);
+        CHECK_ULP_CLOSE(tlo*tlo/2, qh_aos(tlo), 16);
+        CHECK_ULP_CLOSE(thi*thi/2, qh_aos(thi), 16);
+        CHECK_ULP_CLOSE(tlo, qh.prime(tlo), 16);
+        CHECK_ULP_CLOSE(thi, qh.prime(thi), 64);
+        CHECK_ULP_CLOSE(tlo, qh_aos.prime(tlo), 16);
+        CHECK_ULP_CLOSE(thi, qh_aos.prime(thi), 64);
+
+        tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
+        thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
+    }
 }
 
 template<typename Real>