diff --git a/nucleons.jl b/nucleons.jl
index 58d5ecc..32ce4a8 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -62,7 +62,7 @@ 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, Φ0, W0, B0, A0, r_max; dtype=Float64, algo=RK4())
+function boundaryValueFunc(κ, p, Φ0, W0, B0, A0, r_max; dtype=Float64, algo=Tsit5())
     prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, r_max))
     func(E) = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[r_max], save_idxs=[1])[1, 1]
     return func