diff --git a/Project.toml b/Project.toml
index 47f0b481a..b94898885 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,7 +13,6 @@ EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 FastPower = "a4df4552-cc26-4903-aec0-212e50a0e84b"
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -21,7 +20,6 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
-PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
@@ -40,6 +38,7 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
 GTPSA = "b27dd330-f138-47c5-815b-40db9dd9b6e8"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
@@ -55,6 +54,7 @@ DiffEqBaseCUDAExt = "CUDA"
 DiffEqBaseChainRulesCoreExt = "ChainRulesCore"
 DiffEqBaseDistributionsExt = "Distributions"
 DiffEqBaseEnzymeExt = ["ChainRulesCore", "Enzyme"]
+DiffEqBaseForwardDiffExt = ["ForwardDiff"]
 DiffEqBaseGeneralizedGeneratedExt = "GeneralizedGenerated"
 DiffEqBaseGTPSAExt = "GTPSA"
 DiffEqBaseMPIExt = "MPI"
@@ -92,7 +92,6 @@ Measurements = "2"
 MonteCarloMeasurements = "1"
 MuladdMacro = "0.2.1"
 Parameters = "0.12.0"
-PreallocationTools = "0.4"
 PrecompileTools = "1"
 Printf = "1.9"
 RecursiveArrayTools = "3"
diff --git a/src/forwarddiff.jl b/ext/DiffEqBaseForwardDiffExt.jl
similarity index 72%
rename from src/forwarddiff.jl
rename to ext/DiffEqBaseForwardDiffExt.jl
index ec3612fa3..f3bc40a7f 100644
--- a/src/forwarddiff.jl
+++ b/ext/DiffEqBaseForwardDiffExt.jl
@@ -1,15 +1,115 @@
-const DUALCHECK_RECURSION_MAX = 10
+module DiffEqBaseForwardDiffExt
+
+using DiffEqBase, ForwardDiff
+using DiffEqBase.ArrayInterface
+using DiffEqBase: Void, FunctionWrappersWrappers, OrdinaryDiffEqTag, AbstractTimeseriesSolution,
+    RecursiveArrayTools, reduce_tup, _promote_tspan, has_continuous_callback
+import DiffEqBase: hasdualpromote, wrapfun_oop, wrapfun_iip, prob2dtmin,
+                   promote_tspan, anyeltypedual, isdualtype, value, ODE_DEFAULT_NORM,
+                   InternalITP, nextfloat_tdir, DualEltypeChecker, sse
+
+eltypedual(x) = eltype(x) <: ForwardDiff.Dual
+isdualtype(::Type{<:ForwardDiff.Dual}) = true
+const dualT = ForwardDiff.Dual{ForwardDiff.Tag{OrdinaryDiffEqTag, Float64}, Float64, 1}
+dualgen(::Type{T}) where {T} = ForwardDiff.Dual{ForwardDiff.Tag{OrdinaryDiffEqTag, T}, T, 1}
+
+# Copy of the other prob2dtmin dispatch, just for optionality
+function prob2dtmin(tspan, ::ForwardDiff.Dual, use_end_time)
+    t1, t2 = tspan
+    isfinite(t1) || throw(ArgumentError("t0 in the tspan `(t0, t1)` must be finite"))
+    if use_end_time && isfinite(t2 - t1)
+        return max(eps(t2), eps(t1))
+    else
+        return max(eps(typeof(t1)), eps(t1))
+    end
+end
 
-"""
-    promote_dual(::Type{T},::Type{T2})
+function hasdualpromote(u0, t::Number)
+    hasmethod(ArrayInterface.promote_eltype,
+            Tuple{Type{typeof(u0)}, Type{dualgen(eltype(u0))}}) &&
+        hasmethod(promote_rule,
+            Tuple{Type{eltype(u0)}, Type{dualgen(eltype(u0))}}) &&
+        hasmethod(promote_rule,
+            Tuple{Type{eltype(u0)}, Type{typeof(t)}})
+end
+
+const NORECOMPILE_IIP_SUPPORTED_ARGS = (
+    Tuple{Vector{Float64}, Vector{Float64},
+        Vector{Float64}, Float64},
+    Tuple{Vector{Float64}, Vector{Float64},
+        SciMLBase.NullParameters, Float64})
+
+const oop_arglists = (Tuple{Vector{Float64}, Vector{Float64}, Float64},
+    Tuple{Vector{Float64}, SciMLBase.NullParameters, Float64},
+    Tuple{Vector{Float64}, Vector{Float64}, dualT},
+    Tuple{Vector{dualT}, Vector{Float64}, Float64},
+    Tuple{Vector{dualT}, SciMLBase.NullParameters, Float64},
+    Tuple{Vector{Float64}, SciMLBase.NullParameters, dualT})
+
+const NORECOMPILE_OOP_SUPPORTED_ARGS = (Tuple{Vector{Float64},
+        Vector{Float64}, Float64},
+    Tuple{Vector{Float64},
+        SciMLBase.NullParameters, Float64})
+const oop_returnlists = (Vector{Float64}, Vector{Float64},
+    ntuple(x -> Vector{dualT}, length(oop_arglists) - 2)...)
+
+function wrapfun_oop(ff, inputs::Tuple = ())
+    if !isempty(inputs)
+        IT = Tuple{map(typeof, inputs)...}
+        if IT ∉ NORECOMPILE_OOP_SUPPORTED_ARGS
+            throw(NoRecompileArgumentError(IT))
+        end
+    end
+    FunctionWrappersWrappers.FunctionWrappersWrapper(ff, oop_arglists,
+        oop_returnlists)
+end
 
+function wrapfun_iip(ff,
+        inputs::Tuple{T1, T2, T3, T4}) where {T1, T2, T3, T4}
+    T = eltype(T2)
+    dualT = dualgen(T)
+    dualT1 = ArrayInterface.promote_eltype(T1, dualT)
+    dualT2 = ArrayInterface.promote_eltype(T2, dualT)
+    dualT4 = dualgen(promote_type(T, T4))
+
+    iip_arglists = (Tuple{T1, T2, T3, T4},
+        Tuple{dualT1, dualT2, T3, T4},
+        Tuple{dualT1, T2, T3, dualT4},
+        Tuple{dualT1, dualT2, T3, dualT4})
+
+    iip_returnlists = ntuple(x -> Nothing, 4)
+
+    fwt = map(iip_arglists, iip_returnlists) do A, R
+        FunctionWrappersWrappers.FunctionWrappers.FunctionWrapper{R, A}(Void(ff))
+    end
+    FunctionWrappersWrappers.FunctionWrappersWrapper{typeof(fwt), false}(fwt)
+end
+
+const iip_arglists_default = (
+    Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64},
+        Float64},
+    Tuple{Vector{Float64}, Vector{Float64},
+        SciMLBase.NullParameters,
+        Float64
+    },
+    Tuple{Vector{dualT}, Vector{Float64}, Vector{Float64}, dualT},
+    Tuple{Vector{dualT}, Vector{dualT}, Vector{Float64}, dualT},
+    Tuple{Vector{dualT}, Vector{dualT}, Vector{Float64}, Float64},
+    Tuple{Vector{dualT}, Vector{dualT}, SciMLBase.NullParameters,
+        Float64
+    },
+    Tuple{Vector{dualT}, Vector{Float64},
+        SciMLBase.NullParameters, dualT
+    })
+const iip_returnlists_default = ntuple(x -> Nothing, length(iip_arglists_default))
+
+function wrapfun_iip(@nospecialize(ff))
+    fwt = map(iip_arglists_default, iip_returnlists_default) do A, R
+        FunctionWrappersWrappers.FunctionWrappers.FunctionWrapper{R, A}(Void(ff))
+    end
+    FunctionWrappersWrappers.FunctionWrappersWrapper{typeof(fwt), false}(fwt)
+end
 
-Is like the number promotion system, but always prefers a dual number type above
-anything else. For higher order differentiation, it returns the most dualiest of
-them all. This is then used to promote `u0` into the suspected highest differentiation
-space for solving the equation.
-"""
-promote_dual(::Type{T}, ::Type{T2}) where {T, T2} = T
 promote_dual(::Type{T}, ::Type{T2}) where {T <: ForwardDiff.Dual, T2} = T
 function promote_dual(::Type{T},
         ::Type{T2}) where {T <: ForwardDiff.Dual, T2 <: ForwardDiff.Dual}
@@ -38,6 +138,17 @@ function promote_dual(::Type{T},
     ForwardDiff.Dual{T3, promote_dual(V, V2), N}
 end
 
+"""
+    promote_dual(::Type{T},::Type{T2})
+
+
+Is like the number promotion system, but always prefers a dual number type above
+anything else. For higher order differentiation, it returns the most dualiest of
+them all. This is then used to promote `u0` into the suspected highest differentiation
+space for solving the equation.
+"""
+promote_dual(::Type{T}, ::Type{T2}) where {T, T2} = T
+
 # `reduce` and `map` are specialized on tuples to be unrolled (via recursion)
 # Therefore, they can be type stable even with heterogeneous input types.
 # We also don't care about allocating any temporaries with them, as it should
@@ -54,15 +165,12 @@ end
 # For other container types, we probably just want to call `mapreduce`
 @inline diffeqmapreduce(f::F, op::OP, x) where {F, OP} = mapreduce(f, op, x, init = Any)
 
-struct DualEltypeChecker{T, T2}
-    x::T
-    counter::T2
-end
-
 getval(::Val{I}) where {I} = I
 getval(::Type{Val{I}}) where {I} = I
 getval(I::Int) = I
 
+const DUALCHECK_RECURSION_MAX = 10
+
 function (dec::DualEltypeChecker)(::Val{Y}) where {Y}
     isdefined(dec.x, Y) || return Any
     getval(dec.counter) >= DUALCHECK_RECURSION_MAX && return Any
@@ -92,9 +200,9 @@ upconversion is not done automatically, the user is required to upconvert all in
 themselves, for an example of how this can be confusing to a user see
 https://discourse.julialang.org/t/typeerror-in-julia-turing-when-sampling-for-a-forced-differential-equation/82937
 """
-@generated function anyeltypedual(x, ::Type{Val{counter}} = Val{0}) where {counter}
+@generated function anyeltypedual(x, ::Type{Val{counter}}) where {counter}
     x = x.name === Core.Compiler.typename(Type) ? x.parameters[1] : x
-    if x <: ForwardDiff.Dual
+    if isdualtype(x)
         :($x)
     elseif fieldnames(x) === ()
         :(Any)
@@ -180,17 +288,6 @@ function anyeltypedual(::Type{Union{}})
     throw(ForwardDiffAutomaticDetectionFailure())
 end
 
-# Opt out since these are using for preallocation, not differentiation
-function anyeltypedual(x::Union{ForwardDiff.AbstractConfig, Module},
-        ::Type{Val{counter}} = Val{0}) where {counter}
-    Any
-end
-function anyeltypedual(x::Type{T},
-        ::Type{Val{counter}} = Val{0}) where {counter} where {T <:
-                                                              ForwardDiff.AbstractConfig}
-    Any
-end
-
 function anyeltypedual(::Type{<:AbstractTimeseriesSolution{T, N}},
         ::Type{Val{counter}} = Val{0}) where {T, N, counter}
     anyeltypedual(T)
@@ -212,15 +309,6 @@ function anyeltypedual(
     return anyeltypedual((uType, pType), Val{counter})
 end
 
-function anyeltypedual(x::ForwardDiff.DiffResults.DiffResult,
-        ::Type{Val{counter}} = Val{0}) where {counter}
-    Any
-end
-function anyeltypedual(x::Type{T},
-        ::Type{Val{counter}} = Val{0}) where {counter} where {T <:
-                                                              ForwardDiff.DiffResults.DiffResult}
-    Any
-end
 function anyeltypedual(x::SciMLBase.RecipesBase.AbstractPlot,
         ::Type{Val{counter}} = Val{0}) where {counter}
     Any
@@ -229,13 +317,6 @@ function anyeltypedual(x::Returns, ::Type{Val{counter}} = Val{0}) where {counter
     anyeltypedual(x.value, Val{counter})
 end
 
-if isdefined(PreallocationTools, :FixedSizeDiffCache)
-    function anyeltypedual(x::PreallocationTools.FixedSizeDiffCache,
-            ::Type{Val{counter}} = Val{0}) where {counter}
-        Any
-    end
-end
-
 Base.@assume_effects :foldable function __anyeltypedual(::Type{T}) where {T}
     if T isa Union
         promote_dual(anyeltypedual(T.a), anyeltypedual(T.b))
@@ -248,10 +329,7 @@ end
 function anyeltypedual(::Type{T}, ::Type{Val{counter}} = Val{0}) where {counter} where {T}
     __anyeltypedual(T)
 end
-function anyeltypedual(::Type{T},
-        ::Type{Val{counter}} = Val{0}) where {counter} where {T <: ForwardDiff.Dual}
-    T
-end
+
 function anyeltypedual(::Type{T},
         ::Type{Val{counter}} = Val{0}) where {counter} where {T <:
                                                               Union{AbstractArray, Set}}
@@ -371,65 +449,30 @@ end
 anyeltypedual(::@Kwargs{}, ::Type{Val{counter}} = Val{0}) where {counter} = Any
 anyeltypedual(::Type{@Kwargs{}}, ::Type{Val{counter}} = Val{0}) where {counter} = Any
 
-@inline promote_u0(::Nothing, p, t0) = nothing
-
-@inline function promote_u0(u0, p, t0)
-    if SciMLStructures.isscimlstructure(p)
-        _p = SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]
-        if !isequal(_p, p)
-            return promote_u0(u0, _p, t0)
-        end
-    end
-    Tu = eltype(u0)
-    if Tu <: ForwardDiff.Dual
-        return u0
-    end
-    Tp = anyeltypedual(p)
-    if Tp == Any
-        Tp = Tu
-    end
-    Tt = anyeltypedual(t0)
-    if Tt == Any
-        Tt = Tu
-    end
-    Tcommon = promote_type(Tu, Tp, Tt)
-    return if Tcommon <: ForwardDiff.Dual
-        Tcommon.(u0)
-    else
-        u0
-    end
+# Opt out since these are using for preallocation, not differentiation
+function anyeltypedual(x::Union{ForwardDiff.AbstractConfig, Module},
+        ::Type{Val{counter}} = Val{0}) where {counter}
+    Any
+end
+function anyeltypedual(x::Type{T},
+        ::Type{Val{counter}} = Val{0}) where {counter} where {T <:
+                                                              ForwardDiff.AbstractConfig}
+    Any
 end
 
-@inline function promote_u0(u0::AbstractArray{<:Complex}, p, t0)
-    if SciMLStructures.isscimlstructure(p)
-        _p = SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]
-        if !isequal(_p, p)
-            return promote_u0(u0, _p, t0)
-        end
-    end
-    Tu = real(eltype(u0))
-    if Tu <: ForwardDiff.Dual
-        return u0
-    end
-    Tp = anyeltypedual(p)
-    if Tp == Any
-        Tp = Tu
-    end
-    Tt = anyeltypedual(t0)
-    if Tt == Any
-        Tt = Tu
-    end
-    Tcommon = promote_type(eltype(u0), Tp, Tt)
-    return if real(Tcommon) <: ForwardDiff.Dual
-        Tcommon.(u0)
-    else
-        u0
-    end
+function anyeltypedual(x::ForwardDiff.DiffResults.DiffResult,
+        ::Type{Val{counter}} = Val{0}) where {counter}
+    Any
+end
+function anyeltypedual(x::Type{T},
+        ::Type{Val{counter}} = Val{0}) where {counter} where {T <:
+                                                              ForwardDiff.DiffResults.DiffResult}
+    Any
 end
 
-function promote_tspan(u0::AbstractArray{<:ForwardDiff.Dual}, p,
-        tspan::Tuple{<:ForwardDiff.Dual, <:ForwardDiff.Dual}, prob, kwargs)
-    return _promote_tspan(tspan, kwargs)
+function anyeltypedual(::Type{T},
+        ::Type{Val{counter}} = Val{0}) where {counter} where {T <: ForwardDiff.Dual}
+    T
 end
 
 function promote_tspan(u0::AbstractArray{<:ForwardDiff.Dual}, p, tspan, prob, kwargs)
@@ -446,10 +489,17 @@ function promote_tspan(u0::AbstractArray{<:Complex{<:ForwardDiff.Dual}}, p, tspa
     return _promote_tspan(real(eltype(u0)).(tspan), kwargs)
 end
 
+function promote_tspan(u0::AbstractArray{<:ForwardDiff.Dual}, p,
+        tspan::Tuple{<:ForwardDiff.Dual, <:ForwardDiff.Dual}, prob, kwargs)
+    return _promote_tspan(tspan, kwargs)
+end
+
 value(x::Type{ForwardDiff.Dual{T, V, N}}) where {T, V, N} = V
 value(x::ForwardDiff.Dual) = value(ForwardDiff.value(x))
 
-sse(x::Number) = abs2(x)
+unitfulvalue(x::Type{ForwardDiff.Dual{T, V, N}}) where {T, V, N} = V
+unitfulvalue(x::ForwardDiff.Dual) = unitfulvalue(ForwardDiff.unitfulvalue(x))
+
 sse(x::ForwardDiff.Dual) = sse(ForwardDiff.value(x)) + sum(sse, ForwardDiff.partials(x))
 totallength(x::Number) = 1
 function totallength(x::ForwardDiff.Dual)
@@ -482,6 +532,70 @@ end
 
 # Static Arrays don't support the `init` keyword argument for `sum`
 @inline __sum(f::F, args...; init, kwargs...) where {F} = sum(f, args...; init, kwargs...)
-@inline function __sum(f::F, a::StaticArraysCore.StaticArray...; init, kwargs...) where {F}
+@inline function __sum(
+        f::F, a::DiffEqBase.StaticArraysCore.StaticArray...; init, kwargs...) where {F}
     return mapreduce(f, +, a...; init, kwargs...)
 end
+
+# Differentiation of internal solver
+
+function scalar_nlsolve_ad(prob, alg::InternalITP, args...; kwargs...)
+    f = prob.f
+    p = value(prob.p)
+
+    if prob isa IntervalNonlinearProblem
+        tspan = value(prob.tspan)
+        newprob = IntervalNonlinearProblem(f, tspan, p; prob.kwargs...)
+    else
+        u0 = value(prob.u0)
+        newprob = NonlinearProblem(f, u0, p; prob.kwargs...)
+    end
+
+    sol = solve(newprob, alg, args...; kwargs...)
+
+    uu = sol.u
+    if p isa Number
+        f_p = ForwardDiff.derivative(Base.Fix1(f, uu), p)
+    else
+        f_p = ForwardDiff.gradient(Base.Fix1(f, uu), p)
+    end
+
+    f_x = ForwardDiff.derivative(Base.Fix2(f, p), uu)
+    pp = prob.p
+    sumfun = let f_x′ = -f_x
+        ((fp, p),) -> (fp / f_x′) * ForwardDiff.partials(p)
+    end
+    partials = sum(sumfun, zip(f_p, pp))
+    return sol, partials
+end
+
+function SciMLBase.solve(
+        prob::IntervalNonlinearProblem{uType, iip,
+            <:ForwardDiff.Dual{T, V, P}},
+        alg::InternalITP, args...;
+        kwargs...) where {uType, iip, T, V, P}
+    sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
+    return SciMLBase.build_solution(prob, alg, ForwardDiff.Dual{T, V, P}(sol.u, partials),
+        sol.resid; retcode = sol.retcode,
+        left = ForwardDiff.Dual{T, V, P}(sol.left, partials),
+        right = ForwardDiff.Dual{T, V, P}(sol.right, partials))
+end
+
+function SciMLBase.solve(
+        prob::IntervalNonlinearProblem{uType, iip,
+            <:AbstractArray{
+                <:ForwardDiff.Dual{T,
+                V,
+                P},
+            }},
+        alg::InternalITP, args...;
+        kwargs...) where {uType, iip, T, V, P}
+    sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
+
+    return SciMLBase.build_solution(prob, alg, ForwardDiff.Dual{T, V, P}(sol.u, partials),
+        sol.resid; retcode = sol.retcode,
+        left = ForwardDiff.Dual{T, V, P}(sol.left, partials),
+        right = ForwardDiff.Dual{T, V, P}(sol.right, partials))
+end
+
+end
diff --git a/ext/DiffEqBaseGTPSAExt.jl b/ext/DiffEqBaseGTPSAExt.jl
index 7a27569c0..8214a48e9 100644
--- a/ext/DiffEqBaseGTPSAExt.jl
+++ b/ext/DiffEqBaseGTPSAExt.jl
@@ -13,6 +13,9 @@ end
 value(x::TPS) = scalar(x)
 value(::Type{<:TPS{T}}) where {T} = T
 
+unitfulvalue(x::TPS) = scalar(x)
+unitfulvalue(::Type{<:TPS{T}}) where {T} = T
+
 ODE_DEFAULT_NORM(u::TPS, t) = normTPS(u)
 ODE_DEFAULT_NORM(f::F, u::TPS, t) where {F} = normTPS(f(u))
 
diff --git a/ext/DiffEqBaseMeasurementsExt.jl b/ext/DiffEqBaseMeasurementsExt.jl
index c08c640fd..804cf0036 100644
--- a/ext/DiffEqBaseMeasurementsExt.jl
+++ b/ext/DiffEqBaseMeasurementsExt.jl
@@ -19,6 +19,9 @@ DiffEqBase.promote_u0(u0, p::AbstractArray{<:Measurements.Measurement}, t0) = el
 value(x::Type{Measurements.Measurement{T}}) where {T} = T
 value(x::Measurements.Measurement) = Measurements.value(x)
 
+unitfulvalue(x::Type{Measurements.Measurement{T}}) where {T} = T
+unitfulvalue(x::Measurements.Measurement) = Measurements.value(x)
+
 # Support adaptive steps should be errorless
 @inline function DiffEqBase.ODE_DEFAULT_NORM(
         u::AbstractArray{
diff --git a/ext/DiffEqBaseMonteCarloMeasurementsExt.jl b/ext/DiffEqBaseMonteCarloMeasurementsExt.jl
index aba4c894b..2ddb724d8 100644
--- a/ext/DiffEqBaseMonteCarloMeasurementsExt.jl
+++ b/ext/DiffEqBaseMonteCarloMeasurementsExt.jl
@@ -26,6 +26,8 @@ end
 
 DiffEqBase.value(x::Type{MonteCarloMeasurements.AbstractParticles{T, N}}) where {T, N} = T
 DiffEqBase.value(x::MonteCarloMeasurements.AbstractParticles) = mean(x.particles)
+DiffEqBase.unitfulvalue(x::Type{MonteCarloMeasurements.AbstractParticles{T, N}}) where {T, N} = T
+DiffEqBase.unitfulvalue(x::MonteCarloMeasurements.AbstractParticles) = mean(x.particles)
 
 # Support adaptive steps should be errorless
 @inline function DiffEqBase.ODE_DEFAULT_NORM(
diff --git a/ext/DiffEqBaseReverseDiffExt.jl b/ext/DiffEqBaseReverseDiffExt.jl
index 3919386c7..ed2be1ba9 100644
--- a/ext/DiffEqBaseReverseDiffExt.jl
+++ b/ext/DiffEqBaseReverseDiffExt.jl
@@ -23,6 +23,18 @@ end
 DiffEqBase.value(x::ReverseDiff.TrackedReal) = x.value
 DiffEqBase.value(x::ReverseDiff.TrackedArray) = x.value
 
+
+DiffEqBase.unitfulvalue(x::Type{ReverseDiff.TrackedReal{V, D, O}}) where {V, D, O} = V
+function DiffEqBase.unitfulvalue(x::Type{
+        ReverseDiff.TrackedArray{V, D, N, VA, DA},
+}) where {V, D,
+        N, VA,
+        DA}
+    Array{V, N}
+end
+DiffEqBase.unitfulvalue(x::ReverseDiff.TrackedReal) = x.value
+DiffEqBase.unitfulvalue(x::ReverseDiff.TrackedArray) = x.value
+
 # Force TrackedArray from TrackedReal when reshaping W\b
 DiffEqBase._reshape(v::AbstractVector{<:ReverseDiff.TrackedReal}, siz) = reduce(vcat, v)
 
@@ -41,7 +53,7 @@ function DiffEqBase.promote_u0(u0::AbstractArray{<:ReverseDiff.TrackedReal},
 end
 DiffEqBase.promote_u0(u0, p::ReverseDiff.TrackedArray, t0) = ReverseDiff.track(u0)
 function DiffEqBase.promote_u0(
-        u0, p::ReverseDiff.TrackedArray{T}, t0) where {T <: ForwardDiff.Dual}
+        u0, p::ReverseDiff.TrackedArray{T}, t0) where {T <: ReverseDiff.ForwardDiff.Dual}
     ReverseDiff.track(T.(u0))
 end
 DiffEqBase.promote_u0(u0, p::AbstractArray{<:ReverseDiff.TrackedReal}, t0) = eltype(p).(u0)
diff --git a/ext/DiffEqBaseTrackerExt.jl b/ext/DiffEqBaseTrackerExt.jl
index d9e80521f..406927ded 100644
--- a/ext/DiffEqBaseTrackerExt.jl
+++ b/ext/DiffEqBaseTrackerExt.jl
@@ -15,6 +15,11 @@ DiffEqBase.value(x::Type{Tracker.TrackedArray{T, N, A}}) where {T, N, A} = Array
 DiffEqBase.value(x::Tracker.TrackedReal) = x.data
 DiffEqBase.value(x::Tracker.TrackedArray) = x.data
 
+DiffEqBase.unitfulvalue(x::Type{Tracker.TrackedReal{T}}) where {T} = T
+DiffEqBase.unitfulvalue(x::Type{Tracker.TrackedArray{T, N, A}}) where {T, N, A} = Array{T, N}
+DiffEqBase.unitfulvalue(x::Tracker.TrackedReal) = x.data
+DiffEqBase.unitfulvalue(x::Tracker.TrackedArray) = x.data
+
 DiffEqBase.promote_u0(u0::Tracker.TrackedArray, p::Tracker.TrackedArray, t0) = u0
 function DiffEqBase.promote_u0(u0::AbstractArray{<:Tracker.TrackedReal},
         p::Tracker.TrackedArray, t0)
diff --git a/ext/DiffEqBaseUnitfulExt.jl b/ext/DiffEqBaseUnitfulExt.jl
index d74525804..f17c63f65 100644
--- a/ext/DiffEqBaseUnitfulExt.jl
+++ b/ext/DiffEqBaseUnitfulExt.jl
@@ -12,6 +12,10 @@ end
 # Support adaptive errors should be errorless for exponentiation
 value(x::Type{Unitful.AbstractQuantity{T, D, U}}) where {T, D, U} = T
 value(x::Unitful.AbstractQuantity) = x.val
+
+unitfulvalue(x::Type{T}) where {T <: Unitful.AbstractQuantity} = T
+unitfulvalue(x::Unitful.AbstractQuantity) = x
+
 @inline function DiffEqBase.ODE_DEFAULT_NORM(
         u::AbstractArray{
             <:Unitful.AbstractQuantity,
diff --git a/src/DiffEqBase.jl b/src/DiffEqBase.jl
index 49bd50256..9e5c15598 100644
--- a/src/DiffEqBase.jl
+++ b/src/DiffEqBase.jl
@@ -7,7 +7,7 @@ end
 import PrecompileTools
 
 import FastPower
-@deprecate fastpow(x,y) FastPower.fastpower(x,y)
+@deprecate fastpow(x, y) FastPower.fastpower(x, y)
 
 using ArrayInterface
 
@@ -32,8 +32,6 @@ import TruncatedStacktraces
 
 using Setfield
 
-using ForwardDiff
-
 using EnumX
 
 using Markdown
@@ -41,9 +39,6 @@ using Markdown
 using ConcreteStructs: @concrete
 using FastClosures: @closure
 
-# Could be made optional/glue
-import PreallocationTools
-
 import FunctionWrappersWrappers
 
 using SciMLBase
@@ -111,6 +106,14 @@ Reexport.@reexport using SciMLBase
 
 SciMLBase.isfunctionwrapper(x::FunctionWrapper) = true
 
+## Extension Functions
+
+eltypedual(x) = false
+promote_u0(::Nothing, p, t0) = nothing
+isdualtype(::Type{T}) where {T} = false
+
+## Types
+
 """
 $(TYPEDEF)
 """
@@ -132,14 +135,12 @@ include("utils.jl")
 include("stats.jl")
 include("calculate_residuals.jl")
 include("tableaus.jl")
-include("internal_falsi.jl")
 include("internal_itp.jl")
 
 include("callbacks.jl")
 include("common_defaults.jl")
 include("solve.jl")
 include("internal_euler.jl")
-include("forwarddiff.jl")
 include("termination_conditions_deprecated.jl")  # TODO: remove in the next major release
 include("termination_conditions.jl")
 include("norecompile.jl")
diff --git a/src/internal_falsi.jl b/src/internal_falsi.jl
deleted file mode 100644
index f71978bd9..000000000
--- a/src/internal_falsi.jl
+++ /dev/null
@@ -1,158 +0,0 @@
-"""
-  prevfloat_tdir(x, x0, x1)
-
-Move `x` one floating point towards x0.
-"""
-function prevfloat_tdir(x, x0, x1)
-    x1 > x0 ? prevfloat(x) : nextfloat(x)
-end
-
-function nextfloat_tdir(x, x0, x1)
-    x1 > x0 ? nextfloat(x) : prevfloat(x)
-end
-
-function max_tdir(a, b, x0, x1)
-    x1 > x0 ? max(a, b) : min(a, b)
-end
-
-"""
-`InternalFalsi`: A non-allocating regula falsi method, internal to DiffEqBase for
-simpler dependencies.
-"""
-struct InternalFalsi end
-
-function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::InternalFalsi, args...;
-        maxiters = 1000,
-        kwargs...)
-    f = Base.Fix2(prob.f, prob.p)
-    left, right = prob.tspan
-    fl, fr = f(left), f(right)
-
-    if iszero(fl)
-        return SciMLBase.build_solution(prob, alg, left, fl;
-            retcode = ReturnCode.ExactSolutionLeft, left = left,
-            right = right)
-    end
-
-    if iszero(fr)
-        return SciMLBase.build_solution(prob, alg, right, fr;
-            retcode = ReturnCode.ExactSolutionLeft, left = left,
-            right = right)
-    end
-
-    i = 1
-    using_falsi_steps = true
-    while i < maxiters
-        # First, perform a regula falsi iteration
-        if using_falsi_steps
-            if nextfloat_tdir(left, prob.tspan...) == right
-                return SciMLBase.build_solution(prob, alg, left, fl;
-                    retcode = ReturnCode.FloatingPointLimit,
-                    left = left, right = right)
-            end
-            mid = (fr * left - fl * right) / (fr - fl)
-            for i in 1:10
-                mid = max_tdir(left, prevfloat_tdir(mid, prob.tspan...), prob.tspan...)
-            end
-            if mid == right || mid == left
-                using_falsi_steps = false
-                continue
-            end
-            fm = f(mid)
-            if iszero(fm)
-                right = mid
-                using_falsi_steps = false
-                continue
-            end
-            if sign(fl) == sign(fm)
-                fl = fm
-                left = mid
-            else
-                fr = fm
-                right = mid
-            end
-            i += 1
-        end
-
-        # Then, perform a bisection iteration
-        mid = (left + right) / 2
-        (mid == left || mid == right) &&
-            return SciMLBase.build_solution(prob, alg, left, fl;
-                retcode = ReturnCode.FloatingPointLimit,
-                left = left, right = right)
-        fm = f(mid)
-        if iszero(fm)
-            right = mid
-            fr = fm
-        elseif sign(fm) == sign(fl)
-            left = mid
-            fl = fm
-        else
-            right = mid
-            fr = fm
-        end
-        i += 1
-    end
-
-    return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters,
-        left = left, right = right)
-end
-
-function scalar_nlsolve_ad(prob, alg::InternalFalsi, args...; kwargs...)
-    f = prob.f
-    p = value(prob.p)
-
-    if prob isa IntervalNonlinearProblem
-        tspan = value(prob.tspan)
-        newprob = IntervalNonlinearProblem(f, tspan, p; prob.kwargs...)
-    else
-        u0 = value(prob.u0)
-        newprob = NonlinearProblem(f, u0, p; prob.kwargs...)
-    end
-
-    sol = solve(newprob, alg, args...; kwargs...)
-
-    uu = sol.u
-    if p isa Number
-        f_p = ForwardDiff.derivative(Base.Fix1(f, uu), p)
-    else
-        f_p = ForwardDiff.gradient(Base.Fix1(f, uu), p)
-    end
-
-    f_x = ForwardDiff.derivative(Base.Fix2(f, p), uu)
-    pp = prob.p
-    sumfun = let f_x′ = -f_x
-        ((fp, p),) -> (fp / f_x′) * ForwardDiff.partials(p)
-    end
-    partials = sum(sumfun, zip(f_p, pp))
-    return sol, partials
-end
-
-function SciMLBase.solve(
-        prob::IntervalNonlinearProblem{uType, iip,
-            <:ForwardDiff.Dual{T, V, P}},
-        alg::InternalFalsi, args...;
-        kwargs...) where {uType, iip, T, V, P}
-    sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
-    return SciMLBase.build_solution(prob, alg, ForwardDiff.Dual{T, V, P}(sol.u, partials),
-        sol.resid; retcode = sol.retcode,
-        left = ForwardDiff.Dual{T, V, P}(sol.left, partials),
-        right = ForwardDiff.Dual{T, V, P}(sol.right, partials))
-end
-
-function SciMLBase.solve(
-        prob::IntervalNonlinearProblem{uType, iip,
-            <:AbstractArray{
-                <:ForwardDiff.Dual{T,
-                V,
-                P},
-            }},
-        alg::InternalFalsi, args...;
-        kwargs...) where {uType, iip, T, V, P}
-    sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
-
-    return SciMLBase.build_solution(prob, alg, ForwardDiff.Dual{T, V, P}(sol.u, partials),
-        sol.resid; retcode = sol.retcode,
-        left = ForwardDiff.Dual{T, V, P}(sol.left, partials),
-        right = ForwardDiff.Dual{T, V, P}(sol.right, partials))
-end
diff --git a/src/internal_itp.jl b/src/internal_itp.jl
index a5776e0fd..ebdeb9292 100644
--- a/src/internal_itp.jl
+++ b/src/internal_itp.jl
@@ -1,3 +1,20 @@
+"""
+  prevfloat_tdir(x, x0, x1)
+
+Move `x` one floating point towards x0.
+"""
+function prevfloat_tdir(x, x0, x1)
+    x1 > x0 ? prevfloat(x) : nextfloat(x)
+end
+
+function nextfloat_tdir(x, x0, x1)
+    x1 > x0 ? nextfloat(x) : prevfloat(x)
+end
+
+function max_tdir(a, b, x0, x1)
+    x1 > x0 ? max(a, b) : min(a, b)
+end
+
 """
 `InternalITP`: A non-allocating ITP method, internal to DiffEqBase for
 simpler dependencies.
@@ -96,62 +113,3 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem{IP, Tuple{T, T}}, alg::I
     return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters,
         left = left, right = right)
 end
-
-function scalar_nlsolve_ad(prob, alg::InternalITP, args...; kwargs...)
-    f = prob.f
-    p = value(prob.p)
-
-    if prob isa IntervalNonlinearProblem
-        tspan = value(prob.tspan)
-        newprob = IntervalNonlinearProblem(f, tspan, p; prob.kwargs...)
-    else
-        u0 = value(prob.u0)
-        newprob = NonlinearProblem(f, u0, p; prob.kwargs...)
-    end
-
-    sol = solve(newprob, alg, args...; kwargs...)
-
-    uu = sol.u
-    if p isa Number
-        f_p = ForwardDiff.derivative(Base.Fix1(f, uu), p)
-    else
-        f_p = ForwardDiff.gradient(Base.Fix1(f, uu), p)
-    end
-
-    f_x = ForwardDiff.derivative(Base.Fix2(f, p), uu)
-    pp = prob.p
-    sumfun = let f_x′ = -f_x
-        ((fp, p),) -> (fp / f_x′) * ForwardDiff.partials(p)
-    end
-    partials = sum(sumfun, zip(f_p, pp))
-    return sol, partials
-end
-
-function SciMLBase.solve(
-        prob::IntervalNonlinearProblem{uType, iip,
-            <:ForwardDiff.Dual{T, V, P}},
-        alg::InternalITP, args...;
-        kwargs...) where {uType, iip, T, V, P}
-    sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
-    return SciMLBase.build_solution(prob, alg, ForwardDiff.Dual{T, V, P}(sol.u, partials),
-        sol.resid; retcode = sol.retcode,
-        left = ForwardDiff.Dual{T, V, P}(sol.left, partials),
-        right = ForwardDiff.Dual{T, V, P}(sol.right, partials))
-end
-
-function SciMLBase.solve(
-        prob::IntervalNonlinearProblem{uType, iip,
-            <:AbstractArray{
-                <:ForwardDiff.Dual{T,
-                V,
-                P},
-            }},
-        alg::InternalITP, args...;
-        kwargs...) where {uType, iip, T, V, P}
-    sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
-
-    return SciMLBase.build_solution(prob, alg, ForwardDiff.Dual{T, V, P}(sol.u, partials),
-        sol.resid; retcode = sol.retcode,
-        left = ForwardDiff.Dual{T, V, P}(sol.left, partials),
-        right = ForwardDiff.Dual{T, V, P}(sol.right, partials))
-end
diff --git a/src/norecompile.jl b/src/norecompile.jl
index a097fa0f0..273fb14ce 100644
--- a/src/norecompile.jl
+++ b/src/norecompile.jl
@@ -1,29 +1,5 @@
 struct OrdinaryDiffEqTag end
 
-const dualT = ForwardDiff.Dual{ForwardDiff.Tag{OrdinaryDiffEqTag, Float64}, Float64, 1}
-
-dualgen(::Type{T}) where {T} = ForwardDiff.Dual{ForwardDiff.Tag{OrdinaryDiffEqTag, T}, T, 1}
-
-const NORECOMPILE_IIP_SUPPORTED_ARGS = (
-    Tuple{Vector{Float64}, Vector{Float64},
-        Vector{Float64}, Float64},
-    Tuple{Vector{Float64}, Vector{Float64},
-        SciMLBase.NullParameters, Float64})
-
-const oop_arglists = (Tuple{Vector{Float64}, Vector{Float64}, Float64},
-    Tuple{Vector{Float64}, SciMLBase.NullParameters, Float64},
-    Tuple{Vector{Float64}, Vector{Float64}, dualT},
-    Tuple{Vector{dualT}, Vector{Float64}, Float64},
-    Tuple{Vector{dualT}, SciMLBase.NullParameters, Float64},
-    Tuple{Vector{Float64}, SciMLBase.NullParameters, dualT})
-
-const NORECOMPILE_OOP_SUPPORTED_ARGS = (Tuple{Vector{Float64},
-        Vector{Float64}, Float64},
-    Tuple{Vector{Float64},
-        SciMLBase.NullParameters, Float64})
-const oop_returnlists = (Vector{Float64}, Vector{Float64},
-    ntuple(x -> Vector{dualT}, length(oop_arglists) - 2)...)
-
 const NORECOMPILE_ARGUMENT_MESSAGE = """
                                      No-recompile mode is only supported for state arguments
                                      of type `Vector{Float64}`, time arguments of `Float64`
@@ -41,63 +17,15 @@ function Base.showerror(io::IO, e::NoRecompileArgumentError)
     print(io, e.args)
 end
 
-function wrapfun_oop(ff, inputs::Tuple = ())
-    if !isempty(inputs)
-        IT = Tuple{map(typeof, inputs)...}
-        if IT ∉ NORECOMPILE_OOP_SUPPORTED_ARGS
-            throw(NoRecompileArgumentError(IT))
-        end
-    end
-    FunctionWrappersWrappers.FunctionWrappersWrapper(ff, oop_arglists,
-        oop_returnlists)
-end
-
-function wrapfun_iip(ff,
-        inputs::Tuple{T1, T2, T3, T4}) where {T1, T2, T3, T4}
-    T = eltype(T2)
-    dualT = dualgen(T)
-    dualT1 = ArrayInterface.promote_eltype(T1, dualT)
-    dualT2 = ArrayInterface.promote_eltype(T2, dualT)
-    dualT4 = dualgen(promote_type(T, T4))
-
-    iip_arglists = (Tuple{T1, T2, T3, T4},
-        Tuple{dualT1, dualT2, T3, T4},
-        Tuple{dualT1, T2, T3, dualT4},
-        Tuple{dualT1, dualT2, T3, dualT4})
-
-    iip_returnlists = ntuple(x -> Nothing, 4)
-
-    fwt = map(iip_arglists, iip_returnlists) do A, R
-        FunctionWrappersWrappers.FunctionWrappers.FunctionWrapper{R, A}(Void(ff))
-    end
-    FunctionWrappersWrappers.FunctionWrappersWrapper{typeof(fwt), false}(fwt)
+function unwrap_fw(fw::FunctionWrapper)
+    fw.obj[]
 end
 
-const iip_arglists_default = (
-    Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64},
-        Float64},
-    Tuple{Vector{Float64}, Vector{Float64},
-        SciMLBase.NullParameters,
-        Float64
-    },
-    Tuple{Vector{dualT}, Vector{Float64}, Vector{Float64}, dualT},
-    Tuple{Vector{dualT}, Vector{dualT}, Vector{Float64}, dualT},
-    Tuple{Vector{dualT}, Vector{dualT}, Vector{Float64}, Float64},
-    Tuple{Vector{dualT}, Vector{dualT}, SciMLBase.NullParameters,
-        Float64
-    },
-    Tuple{Vector{dualT}, Vector{Float64},
-        SciMLBase.NullParameters, dualT
-    })
-const iip_returnlists_default = ntuple(x -> Nothing, length(iip_arglists_default))
-
-function wrapfun_iip(@nospecialize(ff))
-    fwt = map(iip_arglists_default, iip_returnlists_default) do A, R
-        FunctionWrappersWrappers.FunctionWrappers.FunctionWrapper{R, A}(Void(ff))
-    end
-    FunctionWrappersWrappers.FunctionWrappersWrapper{typeof(fwt), false}(fwt)
+# Default dispatch assumes no ForwardDiff, gets added in the new dispatch
+function wrapfun_iip(ff, inputs)
+    FunctionWrappersWrappers.FunctionWrappersWrapper(Void(ff), (typeof(inputs),), (Nothing,))
 end
 
-function unwrap_fw(fw::FunctionWrapper)
-    fw.obj[]
-end
+function wrapfun_oop(ff, inputs)
+    FunctionWrappersWrappers.FunctionWrappersWrapper(ff, (typeof(inputs),), (typeof(inputs[1]),))
+end
\ No newline at end of file
diff --git a/src/solve.jl b/src/solve.jl
index f3eb3e841..787ad68c2 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -574,7 +574,8 @@ function init_up(prob::AbstractDEProblem, sensealg, u0, p, args...; kwargs...)
         if tstops === nothing && has_kwargs(prob)
             tstops = get(prob.kwargs, :tstops, nothing)
         end
-        if !(tstops isa Union{Nothing, AbstractArray, Tuple, Real}) && !SciMLBase.allows_late_binding_tstops(alg)
+        if !(tstops isa Union{Nothing, AbstractArray, Tuple, Real}) &&
+           !SciMLBase.allows_late_binding_tstops(alg)
             throw(LateBindingTstopsNotSupportedError())
         end
         _prob = get_concrete_problem(prob, isadaptive(alg); u0 = u0, p = p, kwargs...)
@@ -1110,7 +1111,8 @@ function solve_up(prob::Union{AbstractDEProblem, NonlinearProblem}, sensealg, u0
         if tstops === nothing && has_kwargs(prob)
             tstops = get(prob.kwargs, :tstops, nothing)
         end
-        if !(tstops isa Union{Nothing, AbstractArray, Tuple, Real}) && !SciMLBase.allows_late_binding_tstops(alg)
+        if !(tstops isa Union{Nothing, AbstractArray, Tuple, Real}) &&
+           !SciMLBase.allows_late_binding_tstops(alg)
             throw(LateBindingTstopsNotSupportedError())
         end
         _prob = get_concrete_problem(prob, isadaptive(alg); u0 = u0, p = p, kwargs...)
@@ -1282,33 +1284,25 @@ function promote_f(f::F, ::Val{specialize}, u0, p, t) where {F, specialize}
         f = @set f.jac_prototype = similar(f.jac_prototype, uElType)
     end
 
-    @static if VERSION >= v"1.8-"
-        f = if f isa ODEFunction && isinplace(f) && !(f.f isa AbstractSciMLOperator) &&
-               # Some reinitialization code still uses NLSolvers stuff which doesn't
-               # properly tag, so opt-out if potentially a mass matrix DAE
-               f.mass_matrix isa UniformScaling &&
-               # Jacobians don't wrap, so just ignore those cases
-               f.jac === nothing &&
-               ((specialize === SciMLBase.AutoSpecialize && eltype(u0) !== Any &&
-                 RecursiveArrayTools.recursive_unitless_eltype(u0) === eltype(u0) &&
-                 one(t) === oneunit(t) &&
-                 hasmethod(ArrayInterface.promote_eltype,
-                     Tuple{Type{typeof(u0)}, Type{dualgen(eltype(u0))}}) &&
-                 hasmethod(promote_rule,
-                     Tuple{Type{eltype(u0)}, Type{dualgen(eltype(u0))}}) &&
-                 hasmethod(promote_rule,
-                     Tuple{Type{eltype(u0)}, Type{typeof(t)}})) ||
-                (specialize === SciMLBase.FunctionWrapperSpecialize &&
-                 !(f.f isa FunctionWrappersWrappers.FunctionWrappersWrapper)))
-            return unwrapped_f(f, wrapfun_iip(f.f, (u0, u0, p, t)))
-        else
-            return f
-        end
+    f = if f isa ODEFunction && isinplace(f) && !(f.f isa AbstractSciMLOperator) &&
+           # Some reinitialization code still uses NLSolvers stuff which doesn't
+           # properly tag, so opt-out if potentially a mass matrix DAE
+           f.mass_matrix isa UniformScaling &&
+           # Jacobians don't wrap, so just ignore those cases
+           f.jac === nothing &&
+           ((specialize === SciMLBase.AutoSpecialize && eltype(u0) !== Any &&
+             RecursiveArrayTools.recursive_unitless_eltype(u0) === eltype(u0) &&
+             one(t) === oneunit(t) && hasdualpromote(u0, t)) ||
+            (specialize === SciMLBase.FunctionWrapperSpecialize &&
+             !(f.f isa FunctionWrappersWrappers.FunctionWrappersWrapper)))
+        return unwrapped_f(f, wrapfun_iip(f.f, (u0, u0, p, t)))
     else
         return f
     end
 end
 
+hasdualpromote(u0, t) = true
+
 function promote_f(f::SplitFunction, ::Val{specialize}, u0, p, t) where {specialize}
     typeof(f.cache) === typeof(u0) && isinplace(f) ? f : remake(f, cache = zero(u0))
 end
@@ -1462,7 +1456,7 @@ function check_prob_alg_pairing(prob, alg)
         throw(ProblemSolverPairingError(prob, alg))
     end
 
-    if isdefined(prob, :u0) && eltype(prob.u0) <: ForwardDiff.Dual &&
+    if isdefined(prob, :u0) && eltypedual(prob.u0) &&
        !SciMLBase.isautodifferentiable(alg)
         throw(DirectAutodiffError())
     end
diff --git a/src/utils.jl b/src/utils.jl
index fa4488e8b..5b94e7206 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,6 +1,8 @@
 # Handled in Extensions
 value(x) = x
+unitfulvalue(x) = x
 isdistribution(u0) = false
+sse(x::Number) = abs2(x)
 
 _vec(v) = vec(v)
 _vec(v::Number) = v
@@ -22,7 +24,7 @@ end
 
 # This functino requires `eps` to exist, which restricts below `<: Real`
 # Example of a failure is Rational
-function prob2dtmin(tspan, ::Union{AbstractFloat, ForwardDiff.Dual}, use_end_time)
+function prob2dtmin(tspan, ::AbstractFloat, use_end_time)
     t1, t2 = tspan
     isfinite(t1) || throw(ArgumentError("t0 in the tspan `(t0, t1)` must be finite"))
     if use_end_time && isfinite(t2 - t1)
@@ -122,3 +124,65 @@ end
     return __apply_termination_internalnorm(Base.Fix1(maximum, abs), u)
 end
 @inline __apply_termination_internalnorm(f::F, u) where {F} = f(u)
+
+struct DualEltypeChecker{T, T2}
+    x::T
+    counter::T2
+end
+
+anyeltypedual(x) = anyeltypedual(x, Val{0})
+anyeltypedual(x, counter) = Any
+
+function promote_u0(u0, p, t0)
+    if SciMLStructures.isscimlstructure(p)
+        _p = SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]
+        if !isequal(_p, p)
+            return promote_u0(u0, _p, t0)
+        end
+    end
+    Tu = eltype(u0)
+    if isdualtype(Tu)
+        return u0
+    end
+    Tp = anyeltypedual(p, Val{0})
+    if Tp == Any
+        Tp = Tu
+    end
+    Tt = anyeltypedual(t0, Val{0})
+    if Tt == Any
+        Tt = Tu
+    end
+    Tcommon = promote_type(Tu, Tp, Tt)
+    return if isdualtype(Tcommon)
+        Tcommon.(u0)
+    else
+        u0
+    end
+end
+
+function promote_u0(u0::AbstractArray{<:Complex}, p, t0)
+    if SciMLStructures.isscimlstructure(p)
+        _p = SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]
+        if !isequal(_p, p)
+            return promote_u0(u0, _p, t0)
+        end
+    end
+    Tu = real(eltype(u0))
+    if isdualtype(Tu)
+        return u0
+    end
+    Tp = anyeltypedual(p, Val{0})
+    if Tp == Any
+        Tp = Tu
+    end
+    Tt = anyeltypedual(t0, Val{0})
+    if Tt == Any
+        Tt = Tu
+    end
+    Tcommon = promote_type(eltype(u0), Tp, Tt)
+    return if isdualtype(real(Tcommon))
+        Tcommon.(u0)
+    else
+        u0
+    end
+end
diff --git a/test/internal_rootfinder.jl b/test/internal_rootfinder.jl
index 323cc25be..a334eb577 100644
--- a/test/internal_rootfinder.jl
+++ b/test/internal_rootfinder.jl
@@ -1,8 +1,8 @@
 using DiffEqBase
-using DiffEqBase: InternalFalsi, InternalITP, IntervalNonlinearProblem
+using DiffEqBase: InternalITP, IntervalNonlinearProblem
 using ForwardDiff
 
-for Rootfinder in (InternalFalsi, InternalITP)
+for Rootfinder in (InternalITP,)
     rf = Rootfinder()
     # From SimpleNonlinearSolve
     f = (u, p) -> u * u - p