diff --git a/nucleons.jl b/nucleons.jl
index 90e5eec..b3fde98 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -1,4 +1,4 @@
-using DifferentialEquations, Interpolations
+using LinearAlgebra, DifferentialEquations, Interpolations
 include("bisection.jl")
 include("common.jl")
 include("system.jl")
@@ -39,11 +39,26 @@ 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
+
+"Initial boundary condition for u(r)=[g(r), f(r)] at r=0"
+init_BC() = [1.0, 1.0] # TODO: Why not [0.0, 1.0]?
+
 "Solve the Dirac equation and return the wave function u(r)=[g(r), f(r)] where
     the solution would be returned as a 2×mesh_size matrix,
     shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max,
     the other parameters are the same from dirac!(...)."
-function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Tsit5())
+function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
     (f1, f2) = optimized_dirac_potentials(p, s)
 
     # partitioning
@@ -60,19 +75,22 @@ function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Tsit5())
     # right partition
     prob = ODEProblem(dirac!, [0, 1], (s.r_mesh.r_max, r_mid))
     sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=right_r)
-    wf_right = reverse(hcat(sol.u...); dims=2)    
-
+    wf_right = reverse(hcat(sol.u...); dims=2)
+    
     # join two segments
-    rescale_factor_g = wf_left[1, end] / wf_right[1, 1]
-    rescale_factor_f = wf_left[2, end] / wf_right[2, 1]
-    @assert isfinite(rescale_factor_g) && isfinite(rescale_factor_f) "Cannot rescale the right partition"
-    isapprox(rescale_factor_g, rescale_factor_f; rtol=0.03) || @warn "Discontinuity between two partitions"
-    wf_right_rescaled = wf_right .* ((rescale_factor_g + rescale_factor_f) / 2)
-    wf = hcat(wf_left[:, 1:(end - 1)], wf_right_rescaled)
+    u1 = wf_left[:, end]
+    u2 = wf_right[:, 1]
+    if norm(u2) < 1e-10
+        @warn "Right partition too small to rescale, setting to zero"
+        wf_right .= 0.0
+    else
+        proj = only(u1' * u2) / norm(u2)^2
+        wf_right .*= proj
+    end
+    wf = hcat(wf_left[:, 1:(end - 1)], wf_right)
 
     if normalize
-        norm = sum(wf .* wf .* transpose(s.r_mesh.w)) # integration by Reimann sum
-        wf = wf ./ sqrt(norm)
+        wf ./= sqrt(sum(wf .* wf .* transpose(s.r_mesh.w)))
     end
 
     return wf
@@ -81,11 +99,11 @@ end
 "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_mesh.r_max/2, algo=Tsit5())
+function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_mesh.r_max/2, algo=Vern9())
     (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_mesh.r_max, r))
+    prob_left  = ODEProblem(dirac!, init_BC(), (0, r))
     function func(E)
+        prob_right = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_mesh.r_max), (s.r_mesh.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]