diff --git a/nucleons.jl b/nucleons.jl
index 7bee35c..0601768 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -69,10 +69,12 @@ function boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max; dtype=Float64, algo=RK4
 end
 
 "Find all bound energies between E_min (=0) and E_max (=mass) where
+    tol_digits determines the precision for root finding and the threshold for identifying duplicates,
     the other parameters are the same from dirac!(...)."
-function findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n))
+function findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5)
     f(E) = boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max)
-    return find_zeros(f, (E_min, E_max))
+    Es = find_zeros(f, (E_min, E_max); xatol=1/10^tol_digits)
+    return unique(E -> round(E; digits=tol_digits), Es)
 end
 
 "Find all orbitals and return two lists containing κ values and corresponding energies for a single species where
diff --git a/test/Pb208_nucleon_basic.jl b/test/Pb208_nucleon_basic.jl
index ae7e7b0..f823055 100644
--- a/test/Pb208_nucleon_basic.jl
+++ b/test/Pb208_nucleon_basic.jl
@@ -23,7 +23,7 @@ E_max = 939
 
 boundEs = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
 
-binding_Es = round.((p ? M_p : M_n) .- boundEs; digits=3) |> unique
+binding_Es = (p ? M_p : M_n) .- boundEs
 println("binding energies = $binding_Es")
 
 Es = collect(E_min:0.5:E_max)