diff --git a/nucleons.jl b/nucleons.jl
index 690d479..402766e 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -87,10 +87,10 @@ function boundaryValueFunc(κ, p, s::system; dtype=Float64, algo=Tsit5())
     return func
 end
 
-"Find all bound energies between E_min (=0) and E_max (=mass) where
+"Find all bound energies between E_min (=850) and E_max (=938) 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, s::system, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5)
+function findEs(κ, p, s::system, E_min=850, E_max=938, tol_digits=5)
     func = boundaryValueFunc(κ, p, s)
     Es = find_zeros(func, (E_min, E_max); xatol=1/10^tol_digits)
     return unique(E -> round(E; digits=tol_digits), Es)
@@ -98,7 +98,7 @@ end
 
 "Find all orbitals and return two lists containing κ values and corresponding energies for a single species where
     the other parameters are defined above"
-function findAllOrbitals(p, s::system, E_min=0, E_max=(p ? M_p : M_n))
+function findAllOrbitals(p, s::system, E_min=850, E_max=938)
     κs = Int[]
     Es = Float64[]
     # start from κ=-1 and go both up and down
@@ -160,7 +160,7 @@ end
 
 "Solve the Dirac equation and calculate scalar and vector densities of a nucleon species where
     the other parameters are defined above"
-function solveNucleonDensity(p, s::system, E_min=800, E_max=939)
+function solveNucleonDensity(p, s::system, E_min=850, E_max=938)
     κs, Es = findAllOrbitals(p, s, E_min, E_max)
     occs = fillNucleons(Z_or_N(s, p), κs, Es)
     return calculateNucleonDensity(κs, Es, occs, p, s)
diff --git a/test/Linear_nucleon_spectrum.jl b/test/Linear_nucleon_spectrum.jl
index 5714e9c..e8084a4 100644
--- a/test/Linear_nucleon_spectrum.jl
+++ b/test/Linear_nucleon_spectrum.jl
@@ -20,10 +20,7 @@ s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-E_min = 860
-E_max = 939
-
-κs, Es = findAllOrbitals(p, s, E_min, E_max)
+κs, Es = findAllOrbitals(p, s)
 Ebinds = (p ? M_p : M_n) .- Es
 occ = fillNucleons(Z_or_N(s, p), κs, Es)
 
diff --git a/test/Pb208_nucleon_basic.jl b/test/Pb208_nucleon_basic.jl
index a5f2edf..9d31bd6 100644
--- a/test/Pb208_nucleon_basic.jl
+++ b/test/Pb208_nucleon_basic.jl
@@ -22,7 +22,7 @@ s.B0 = Rs
 s.A0 = As
 
 E_min = 850
-E_max = 939
+E_max = 938
 
 boundEs = findEs(κ, p, s, E_min, E_max)
 
diff --git a/test/Pb208_nucleon_spectrum.jl b/test/Pb208_nucleon_spectrum.jl
index c12ea6d..0704593 100644
--- a/test/Pb208_nucleon_spectrum.jl
+++ b/test/Pb208_nucleon_spectrum.jl
@@ -22,10 +22,7 @@ s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-E_min = 800
-E_max = 939
-
-κs, Es = findAllOrbitals(p, s, E_min, E_max)
+κs, Es = findAllOrbitals(p, s)
 Ebinds = (p ? M_p : M_n) .- Es
 occ = fillNucleons(Z_or_N(s, p), κs, Es)
 
diff --git a/test/Pb208_nucleon_wf.jl b/test/Pb208_nucleon_wf.jl
index 7898949..0018384 100644
--- a/test/Pb208_nucleon_wf.jl
+++ b/test/Pb208_nucleon_wf.jl
@@ -21,10 +21,7 @@ s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-E_min = 800
-E_max = 939
-
-groundE = findEs(κ, p, s, E_min, E_max) |> minimum
+groundE = findEs(κ, p, s) |> minimum
 println("ground state E = $groundE")
 
 wf = solveNucleonWf(κ, p, groundE, s)