diff --git a/Project.toml b/Project.toml
index 0f0d325..65b716f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,3 +1,4 @@
 [deps]
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
diff --git a/dirac.jl b/dirac.jl
index 024e6ee..6a6a1c8 100644
--- a/dirac.jl
+++ b/dirac.jl
@@ -1,4 +1,4 @@
-using DifferentialEquations
+using DifferentialEquations, Roots
 
 ħc = 197.327 # ħc in MeVfm
 M_n = 939.5654133 # Neutron mass in MeV/c2
@@ -28,4 +28,15 @@ function boundaryValue(κ, M, E, S, V, r_max, r_min=r_max/1000)
     prob = ODEProblem(dirac!, [0, 1], (r_min, r_max))
     sol = solve(prob, RK4(), p=(κ, M, E, S, V))
     return sol(r_max)[1]
-end
\ No newline at end of file
+end
+
+"Find all bound energies between E_min (=0) and E_max (=M) where
+κ is the generalized angular momentum,
+M is the mass in MeV/c2,
+S(r) & V(r) are functions corresponding to scalar and vector potentials in MeV,
+r_max is the outer boundary,
+r_min (=r_max/1000) is inside boundary which cannot be 0 due to the centrifugal term."
+function findEs(κ, M, S, V, r_max, r_min=r_max/1000, E_min=0, E_max=M)
+    f(E) = boundaryValue(κ, M, E, S, V, r_max, r_min)
+    return find_zeros(f, (E_min, E_max))
+end
diff --git a/test/Pb208.jl b/test/Pb208.jl
index a948354..b5368f7 100644
--- a/test/Pb208.jl
+++ b/test/Pb208.jl
@@ -15,8 +15,18 @@ V_interp = linear_interpolation(xs, Vs)
 R_interp = linear_interpolation(xs, Rs)
 A_interp = linear_interpolation(xs, As)
 
-Es = collect(840:0.5:940)
-boundaryVals = [boundaryValue(-1, M_n, E, S_interp, V_interp, maximum(xs))^2 for E in Es]
+κ = -1
+M = M_n
+r_max = maximum(xs)
+E_min = M - 100
+E_max = M
+
+boundEs = findEs(κ, M_n, S_interp, V_interp, r_max, r_max/1000, E_min, E_max)
+println("bound E = $boundEs")
+
+Es = collect(E_min:0.5:E_max)
+boundaryVals = [boundaryValue(κ, M_n, E, S_interp, V_interp, r_max)^2 for E in Es]
 
 plot(Es, boundaryVals, yscale=:log10, label="g(r_max)^2")
+vline!(boundEs, label="bound E")
 xlabel!("E (MeV)")