diff --git a/example/daubechies_wavelets/daubechies_wavelet_plots.cpp b/example/daubechies_wavelets/daubechies_wavelet_plots.cpp
index d4fcb5952..4d898ea4e 100644
--- a/example/daubechies_wavelets/daubechies_wavelet_plots.cpp
+++ b/example/daubechies_wavelets/daubechies_wavelet_plots.cpp
@@ -121,10 +121,10 @@ void plot_condition_number()
     daub.write_all();
 }
 
-template<typename CoarseReal, typename PreciseReal, int p, class PhiPrecise>
-void do_ulp(int coarse_refinements, PhiPrecise phi_precise)
+template<typename CoarseReal, typename PreciseReal, int p, class PsiPrecise>
+void do_ulp(int coarse_refinements, PsiPrecise psi_precise)
 {
-    auto phi_coarse = boost::math::daubechies_wavelet<CoarseReal, p>(coarse_refinements);
+    auto psi_coarse = boost::math::daubechies_wavelet<CoarseReal, p>(coarse_refinements);
 
     std::string title = std::to_string(p) + " vanishing moment ULP plot at " + std::to_string(coarse_refinements) + " refinements and " + boost::core::demangle(typeid(CoarseReal).name()) + " precision";
     title = "";
@@ -134,7 +134,7 @@ void do_ulp(int coarse_refinements, PhiPrecise phi_precise)
     int clip = 20;
     int horizontal_lines = 8;
     int vertical_lines = 2*p - 1;
-    quicksvg::ulp_plot<decltype(phi_coarse), CoarseReal, decltype(phi_precise), PreciseReal>(phi_coarse, phi_precise, CoarseReal(0), phi_coarse.support().second, title, filename, samples, GRAPH_WIDTH, clip, horizontal_lines, vertical_lines);
+    quicksvg::ulp_plot<decltype(psi_coarse), CoarseReal, decltype(psi_precise), PreciseReal>(psi_coarse, psi_precise, CoarseReal(psi_coarse.support().first), psi_coarse.support().second, title, filename, samples, GRAPH_WIDTH, clip, horizontal_lines, vertical_lines);
 }
 
 
@@ -148,7 +148,7 @@ int main()
     using PreciseReal = float128;
     using CoarseReal = double;
     int precise_refinements = 22;
-    constexpr const int p = 8;
+    constexpr const int p = 9;
     std::cout << "Computing precise wavelet function in " << boost::core::demangle(typeid(PreciseReal).name()) << " precision.\n";
     auto phi_precise = boost::math::daubechies_wavelet<PreciseReal, p>(precise_refinements);
     std::cout << "Beginning comparison with functions computed in " << boost::core::demangle(typeid(CoarseReal).name()) << " precision.\n";
diff --git a/example/daubechies_wavelets/wavelet_transform.cpp b/example/daubechies_wavelets/wavelet_transform.cpp
index 8a5c3c959..67cdde803 100644
--- a/example/daubechies_wavelets/wavelet_transform.cpp
+++ b/example/daubechies_wavelets/wavelet_transform.cpp
@@ -7,30 +7,8 @@
 #include <cstdint>
 #include <cmath>
 #include <boost/math/quadrature/wavelet_transforms.hpp>
-#include <boost/gil.hpp>
-#include <boost/gil/extension/io/png.hpp>
-#include <boost/gil/extension/io/jpeg.hpp>
-#include <boost/gil/image.hpp>
 #include <Eigen/Dense>
 
-namespace bg = boost::gil;
-/*
-int main(int argc, char *argv[])
-{
-    int height = 1000;
-    int width = 1.618*height;
-    auto img = bg::gray16_image_t{width, height, bg::gray16_pixel_t {0}};
-    auto view = bg::view(img);
-    auto count = std::uint16_t {0};
-    for (auto it = view.begin(); it != view.end(); ++it)
-    {
-        *it = count++;
-    }
-
-    //bg::write_view("img.png", bg::const_view(img), bg::png_tag());
-    bg::write_view("img.jpeg", bg::const_view(img), bg::png_tag());
-    return 0;
-}*/
 
 int main()
 {
@@ -65,4 +43,4 @@ int main()
     std::cout << "W[f](-s,t) = " << Wf(-s, t) << "\n";
     std::cout << "W[g](-s,t) = " << Wg(-s, t) << "\n";
 
-}
\ No newline at end of file
+}
diff --git a/test/daubechies_wavelet_test.cpp b/test/daubechies_wavelet_test.cpp
index d9b9f0796..f1614d10e 100644
--- a/test/daubechies_wavelet_test.cpp
+++ b/test/daubechies_wavelet_test.cpp
@@ -16,6 +16,7 @@
 #include <boost/hana/ext/std/integer_sequence.hpp>
 #include <boost/math/tools/condition_numbers.hpp>
 #include <boost/math/special_functions/daubechies_wavelet.hpp>
+#include <boost/math/special_functions/next.hpp>
 #include <boost/math/quadrature/trapezoidal.hpp>
 
 #ifdef BOOST_HAS_FLOAT128
@@ -32,17 +33,10 @@ void test_exact_value()
 {
     // The global phase of the wavelet is not constrained by anything other than convention.
     // Make sure that our conventions match the rest of the world:
-
     auto psi = boost::math::daubechies_wavelet<Real, 2>(2);
-    auto phi = boost::math::daubechies_wavelet<Real, 2>(2);
-    auto h = boost::math::filters::daubechies_scaling_filter<Real, 2>();
-
     Real computed = psi(1);
-    // this expression for expected is wrong!
-    Real expected = root_two<Real>()*(-h[0]*phi(1) + h[1]*phi(2));
-    CHECK_ULP_CLOSE(expected, computed, 1);
-
-    std::cout << "psi(" << 1 << ") = " << psi(1) << "\n";
+    Real expected = -1.366025403784439;
+    CHECK_MOLLIFIED_CLOSE(expected, computed, 0.0001);
 }
 
 template<typename Real, int p>
@@ -77,14 +71,33 @@ void test_quadratures()
         std::cerr << "  Error estimate: " << error_estimate << ", L1 norm: " << L1 << "\n";
     }
     // psi is orthogonal to its integer translates: \int \psi(x-k) \psi(x) \, \mathrm{d}x = 0
-    // psi has L2 norm 1:
-
     // g_n = 1/sqrt(2) <psi(t/2), phi(t-n)> (Mallat, 7.55)
+
+    // Now hit the boundary. Much can go wrong here; this just tests for segfaults:
+    int samples = 500;
+    Real xlo = a;
+    Real xhi = b;
+    for (int i = 0; i < samples; ++i)
+    {
+        CHECK_ULP_CLOSE(Real(0), psi(xlo), 0);
+        CHECK_ULP_CLOSE(Real(0), psi(xhi), 0);
+        xlo = std::nextafter(xlo, std::numeric_limits<Real>::lowest());
+        xhi = std::nextafter(xhi, std::numeric_limits<Real>::max());
+    }
+
+    xlo = a;
+    xhi = b;
+    for (int i = 0; i < samples; ++i) {
+        assert(abs(psi(xlo)) <= 5);
+        assert(abs(psi(xhi)) <= 5);
+        xlo = std::nextafter(xlo, std::numeric_limits<Real>::max());
+        xhi = std::nextafter(xhi, std::numeric_limits<Real>::lowest());
+    }
 }
 
 int main()
 {
-    //test_exact_value<double>();
+    test_exact_value<double>();
 
     boost::hana::for_each(std::make_index_sequence<17>(), [&](auto i){
         test_quadratures<float, i+3>();