diff --git a/nucleons.jl b/nucleons.jl
index 271684c..372c717 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -78,13 +78,18 @@ function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true
     return wf
 end
 
-"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 boundaryValueFunc(κ, p::Bool, s::system; dtype=Float64, algo=Tsit5())
-    prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, s.r_max))
+"Returns a function that solves the Dirac equation in two partitions and returns 
+    the determinant of [g_left(r) g_right(r); f_left(r) f_right(r)],
+    where is r is in fm."
+function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Tsit5())
     (f1, f2) = optimized_dirac_potentials(p, s)
-    func(E) = solve(prob, algo, p=(κ, E, f1, f2), saveat=[s.r_max], save_idxs=[1])[1, 1]
+    prob_left  = ODEProblem(dirac!, [0.0, 1.0], (0, r))
+    prob_right = ODEProblem(dirac!, [0.0, 1.0], (s.r_max, r))
+    function func(E)
+        u_left  = solve(prob_left,  algo, p=(κ, E, f1, f2), saveat=[r])
+        u_right = solve(prob_right, algo, p=(κ, E, f1, f2), saveat=[r])
+        return u_left[1, 1] * u_right[2, 1] - u_right[1, 1] * u_left[2, 1]
+    end
     return func
 end
 
@@ -92,7 +97,7 @@ end
     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::Bool, s::system, E_min=850.0, E_max=938.0, tol_digits=5)
-    func = boundaryValueFunc(κ, p, s)
+    func = determinantFunc(κ, p, s)
     Es = find_all_zeros(func, E_min, E_max; partitions=20, tol=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_shooting.jl
similarity index 79%
rename from test/Pb208_nucleon_basic.jl
rename to test/Pb208_shooting.jl
index 1e0de97..05dc7ed 100644
--- a/test/Pb208_nucleon_basic.jl
+++ b/test/Pb208_shooting.jl
@@ -29,10 +29,10 @@ boundEs = findEs(κ, p, s, E_min, E_max)
 binding_Es = (p ? M_p : M_n) .- boundEs
 println("binding energies = $binding_Es")
 
-func = boundaryValueFunc(κ, p, s)
-Es = collect(E_min:0.5:E_max)
-boundaryVals = [func(E)^2 for E in Es]
+func = determinantFunc(κ, p, s)
+Es = collect(E_min:0.1:E_max)
+dets = [func(E)^2 for E in Es]
 
-plot(Es, boundaryVals, yscale=:log10, label="g(r_max)^2")
+plot(Es, dets, yscale=:log10, label="determinant^2")
 vline!(boundEs, label="bound E")
 xlabel!("E (MeV)")