diff --git a/NEWS.md b/NEWS.md
index a7fcf92b4a..50c24df127 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,7 @@
 # ggplot2 (development version)
 
+* The bounded density option in `stat_density()` uses a wider range to
+  prevent discontinuities (#5641).
 * `geom_raster()` now falls back to rendering as `geom_rect()` when coordinates
   are not Cartesian (#5503).
 * `stat_ecdf()` now has an optional `weight` aesthetic (@teunbrand, #5058).
diff --git a/R/stat-density.R b/R/stat-density.R
index 1b84f2cdb4..27cafd4b6d 100644
--- a/R/stat-density.R
+++ b/R/stat-density.R
@@ -148,10 +148,23 @@ compute_density <- function(x, w, from, to, bw = "nrd0", adjust = 1,
   bw <- precompute_bw(x, bw)
   # Decide whether to use boundary correction
   if (any(is.finite(bounds))) {
-    dens <- stats::density(x, weights = w, bw = bw, adjust = adjust,
-                           kernel = kernel, n = n)
+    # To prevent discontinuities, we widen the range before calling the
+    # unbounded estimator (#5641).
+    bounds   <- sort(bounds)
+    range    <- range(from, to)
+    width    <- diff(range)
+    range[1] <- range[1] - width * as.numeric(is.finite(bounds[1]))
+    range[2] <- range[2] + width * as.numeric(is.finite(bounds[2]))
+    n <- n * (sum(is.finite(bounds)) + 1)
 
-    dens <- reflect_density(dens = dens, bounds = bounds, from = from, to = to)
+    dens <- stats::density(
+      x, weights = w, bw = bw, adjust = adjust,
+      kernel = kernel, n = n, from = range[1], to = range[2]
+    )
+    dens <- reflect_density(
+      dens = dens, bounds = bounds,
+      from = range[1], to = range[2]
+    )
   } else {
     dens <- stats::density(x, weights = w, bw = bw, adjust = adjust,
                            kernel = kernel, n = n, from = from, to = to)
diff --git a/tests/testthat/test-stat-density.R b/tests/testthat/test-stat-density.R
index 0274d83d40..8a772b8873 100644
--- a/tests/testthat/test-stat-density.R
+++ b/tests/testthat/test-stat-density.R
@@ -66,7 +66,7 @@ test_that("stat_density uses `bounds`", {
     expect_equal(
       orig_density(test_sample) + left_reflection + right_reflection,
       plot_density(test_sample),
-      tolerance = 1e-4
+      tolerance = 1e-3
     )
   }
 
@@ -94,7 +94,7 @@ test_that("stat_density handles data outside of `bounds`", {
       stat_density(bounds = c(cutoff, Inf))
   )
 
-  expect_equal(data_actual, data_expected)
+  expect_equal(data_actual, data_expected, tolerance = 1e-4)
 })
 
 test_that("compute_density succeeds when variance is zero", {