diff --git a/nucleons.jl b/nucleons.jl
index 63b3d83..fd7b439 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -55,13 +55,13 @@ end
 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
-    divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix,
-    the other parameters are the same from dirac!(...)."
-function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
+    matching_point marks the partitioning position (between 0.0 and 1.0) for shooting method,
+    the solution would be returned as a 2×(1+divs) matrix."
+function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9(), matching_point=0.15)
     (f1, f2) = optimized_dirac_potentials(p, s)
 
     # partitioning
-    mid_idx = s.divs ÷ 2
+    mid_idx = s.divs * matching_point |> round |> Int
     r_mid = rs(s)[mid_idx]
     left_r = rs(s)[1:mid_idx]
     right_r = rs(s)[mid_idx:end]