Use finite-difference AD for BVP default nonlinear solvers

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

Author: ChrisRackauckas

lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl

@@ -6,7 +6,7 @@
 # resolved
 function __FastShortcutBVPCompatibleNLLSPolyalg(
         ::Type{T} = Float64; concrete_jac = nothing,
-        linsolve = nothing, autodiff = nothing, kwargs...
+        linsolve = nothing, autodiff = AutoFiniteDiff(), kwargs...
     ) where {T}
     if T <: Complex
         algs = (
@@ -30,7 +30,7 @@ end
 
 function __FastShortcutBVPCompatibleNonlinearPolyalg(
         ::Type{T} = Float64; concrete_jac = nothing,
-        linsolve = nothing, autodiff = nothing
+        linsolve = nothing, autodiff = AutoFiniteDiff()
     ) where {T}
     if T <: Complex
         algs = (NewtonRaphson(; concrete_jac, linsolve, autodiff),)