diff --git a/nucleons.jl b/nucleons.jl
index b57c66c..58d5ecc 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -59,21 +59,21 @@ function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, n
     return wf
 end
 
-"Solve the Dirac equation and return g(r=r_max) where
+"Returns a function that solves the Dirac equation and returns g(r=r_max) where
     r_max is the outer boundary in fm,
     the other parameters are the same from dirac!(...)."
-function boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max; dtype=Float64, algo=RK4())
+function boundaryValueFunc(κ, p, Φ0, W0, B0, A0, r_max; dtype=Float64, algo=RK4())
     prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, r_max))
-    sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[r_max], save_idxs=[1])
-    return sol[1, 1]
+    func(E) = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[r_max], save_idxs=[1])[1, 1]
+    return func
 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), tol_digits=5)
-    f(E) = boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max)
-    Es = find_zeros(f, (E_min, E_max); xatol=1/10^tol_digits)
+    func = boundaryValueFunc(κ, p, Φ0, W0, B0, A0, r_max)
+    Es = find_zeros(func, (E_min, E_max); xatol=1/10^tol_digits)
     return unique(E -> round(E; digits=tol_digits), Es)
 end
 
diff --git a/test/Pb208_nucleon_basic.jl b/test/Pb208_nucleon_basic.jl
index f823055..481d55b 100644
--- a/test/Pb208_nucleon_basic.jl
+++ b/test/Pb208_nucleon_basic.jl
@@ -26,8 +26,9 @@ boundEs = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_
 binding_Es = (p ? M_p : M_n) .- boundEs
 println("binding energies = $binding_Es")
 
+func = boundaryValueFunc(κ, p, S_interp, V_interp, R_interp, A_interp, r_max)
 Es = collect(E_min:0.5:E_max)
-boundaryVals = [boundaryValue(κ, p, E, S_interp, V_interp, R_interp, A_interp, r_max)^2 for E in Es]
+boundaryVals = [func(E)^2 for E in Es]
 
 plot(Es, boundaryVals, yscale=:log10, label="g(r_max)^2")
 vline!(boundEs, label="bound E")