diff --git a/include/boost/math/interpolators/detail/bezier_polynomial_detail.hpp b/include/boost/math/interpolators/detail/bezier_polynomial_detail.hpp
index 18ce3e20a..60fe5ab5d 100644
--- a/include/boost/math/interpolators/detail/bezier_polynomial_detail.hpp
+++ b/include/boost/math/interpolators/detail/bezier_polynomial_detail.hpp
@@ -73,7 +73,7 @@ public:
             }
         }
 
-        decasteljau_recursion(scratch_space, scratch_space.size(), t);
+        decasteljau_recursion(scratch_space, control_points_.size() - 1, t);
         return scratch_space[0];
     }
 
diff --git a/test/git_issue_959.cpp b/test/git_issue_959.cpp
new file mode 100644
index 000000000..1a6bf48fa
--- /dev/null
+++ b/test/git_issue_959.cpp
@@ -0,0 +1,42 @@
+// Copyright Matt Borland, 2023
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <iostream>
+#include <array>
+#include <vector>
+#include <boost/math/interpolators/bezier_polynomial.hpp>
+
+using tControlPoints = std::vector<std::array<double, 3>>;
+
+void interpolateWithPoints(tControlPoints cp)
+{
+    const auto cpSize = cp.size();
+    auto bp = boost::math::interpolators::bezier_polynomial(std::move(cp));
+
+    // Interpolate at t = 0.5:
+    std::array<double, 3> point = bp(0.5);
+    std::cout << cpSize << " points, t = 0.5:\n";
+
+    for (const auto& c : point) 
+    {
+        std::cout << "  " << c << "\n";
+    }
+}
+
+
+int main(void)
+{
+    auto cp3 = tControlPoints{{0,0,0}, {1,0,0}, {0,1,0}};
+    interpolateWithPoints(cp3);
+
+    auto cp4 = tControlPoints{{0,0,0}, {1,0,0}, {0,1,0}, {0,0,1}};
+    interpolateWithPoints(cp4);
+
+    auto cp3b = tControlPoints{{0,0,0}, {1,0,0}, {0,1,0}};
+    interpolateWithPoints(cp3b);
+
+    return 0;
+}