diff --git a/Project.toml b/Project.toml
index 7f62bd5..0e86277 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
+PolynomialRoots = "3a141323-8675-5d76-9d11-e1df1406c778"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
diff --git a/src/LinearFitting.jl b/src/LinearFitting.jl
index 28e7e94..5584433 100644
--- a/src/LinearFitting.jl
+++ b/src/LinearFitting.jl
@@ -3,35 +3,42 @@ This module provides some linear fitting methods.
 """
 module LinearFitting
 
-using Polynomials: Polynomial, fit, derivative, roots, coeffs
+using Polynomials: Polynomial, fit, derivative, coeffs
+using PolynomialRoots: roots
 
 using ..Collections: PolynomialEOS
 
 export linfit
 
-_islocalminimum(f, x, δx) = f(x) < f(x - δx) && f(x) < f(x + δx)
+_islocalminimum(y, x, δx) = y(x) < y(x - δx) && y(x) < y(x + δx)
 
-function _findglobalminimum(f, localminima, δx)
-    if length(localminima) == 0
-        error("no volume minimizes the energy!")
+function _findlocalminima(y, xs)
+    y′ = derivative(y, 1)
+    δx = minimum(diff(xs)) / 10
+    localminima = eltype(xs)[]
+    for x in real(filter(isreal, roots(coeffs(y′))))  # Complex volumes are meaningless
+        if _islocalminimum(y, x, δx)
+            push!(localminima, x)
+        end
+    end
+    return localminima
+end # function _findlocalminima
+
+function _findglobalminimum(y, localminima)
+    # https://stackoverflow.com/a/21367608/3260253
+    if isempty(localminima)
+        error("no local minima found!")
     else
-        f0, i = findmin(f.(localminima))
+        y0, i = findmin(y.(localminima))
         x0 = localminima[i]
-        return x0, f0
+        return x0, y0
     end
 end # function _findglobalminimum
 
 function linfit(volumes, energies, deg = 3)
     poly = fit(volumes, energies, deg)
-    poly1d = derivative(poly, 1)
-    δx = minimum(diff(volumes)) / 10
-    localminima = eltype(volumes)[]
-    for x in real(filter(isreal, roots(poly1d)))  # Complex volume is meaningless
-        if _islocalminimum(poly, x, δx)
-            push!(localminima, x)
-        end
-    end
-    v0, e0 = _findglobalminimum(poly, localminima, δx)
+    localminima = _findlocalminima(poly, volumes)
+    v0, e0 = _findglobalminimum(poly, localminima)
     bs = Tuple(derivative(poly, n)(v0) / factorial(n) for n in 1:deg)
     return PolynomialEOS(v0, bs, e0)
 end # function linfit