diff --git a/Project.toml b/Project.toml
index af12e4b7a..b2107a678 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,6 +13,7 @@ PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
@@ -48,6 +49,7 @@ PDMats = "0.10, 0.11"
 Printf = "<0.0.1, 1"
 QuadGK = "2"
 Random = "<0.0.1, 1"
+Roots = "2"
 SparseArrays = "<0.0.1, 1"
 SpecialFunctions = "1.2, 2"
 StableRNGs = "1"
diff --git a/src/quantilealgs.jl b/src/quantilealgs.jl
index 8815811b9..bcef0d507 100644
--- a/src/quantilealgs.jl
+++ b/src/quantilealgs.jl
@@ -61,8 +61,7 @@ function newton_impl(Δ, xs::T=mode(d), xrtol::Real=1e-12) where {T}
 end
 
 function newton((f,df), xs::T=mode(d), xrtol::Real=1e-12) where {T}
-    Δ(x) = f(x)/df(x)
-    return newton_impl(Δ, xs, xrtol)
+    return find_zero((f,df), xs, Roots.Newton(), xatol=0, xrtol=xrtol, atol=0, rtol=eps(float(T)), maxiters=typemax(Int))
 end
 
 function quantile_newton(d::ContinuousUnivariateDistribution, p::Real, xs::Real=mode(d), xrtol::Real=1e-12)