diff --git a/nucleons.jl b/nucleons.jl
index 4796bbf..7fbc194 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -55,7 +55,7 @@ end
     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,
     the other parameters are the same from dirac!(...)."
-function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true, algo=Tsit5())
+function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true, algo=Vern9())
     (f1, f2) = optimized_dirac_potentials(p, s)
 
     if shooting
@@ -95,7 +95,7 @@ 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_max/2, algo=Tsit5())
+function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Vern9())
     (f1, f2) = optimized_dirac_potentials(p, s)
     prob_left  = ODEProblem(dirac!, [0.0, 1.0], (0, r))
     function func(E)