diff --git a/nucleons.jl b/nucleons.jl
index 4b92e20..c70a9be 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -39,6 +39,18 @@ function optimized_dirac_potentials(p::Bool, s::system)
     return (f1, f2)
 end
 
+"Approximate boundary condition for u(r)=[g(r), f(r)] at r -> ∞ where
+    κ is the generalized angular momentum,
+    p is true for proton and false for neutron,
+    E is the energy in MeV,
+    r is the radius in fm."
+function asymp_BC(κ::Int, p::Bool, E::Float64, r::Float64)
+    M = p ? M_p : M_n
+    g = 1.0
+    f = ħc / (E + M) * (-√(M^2 - E^2) + κ/r) * g 
+    return [g, f]
+end
+
 "Solve the Dirac equation and return the wave function u(r)=[g(r), f(r)] where
     divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix,
     shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max,
@@ -48,7 +60,7 @@ function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true
 
     if shooting
         @assert s.divs % 2 == 0 "divs must be an even number when shooting=true"
-        prob = ODEProblem(dirac!, [0, 1], (s.r_max, s.r_max / 2))
+        prob = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, s.r_max / 2))
         sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s))
         wf_right = reverse(hcat(sol.u...); dims=2)
         next_r_max = s.r_max / 2 # for the next step
@@ -83,8 +95,8 @@ end
 function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Tsit5())
     (f1, f2) = optimized_dirac_potentials(p, s)
     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)
+        prob_right = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r))
         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]