how about now

Project.toml

@@ -49,7 +49,7 @@ ADTypes = "0.1, 0.2"
 Adapt = "1.1, 2.0, 3.0"
 ArrayInterface = "6, 7"
 DataStructures = "0.18"
-DiffEqBase = "6.128.2"
+DiffEqBase = "6.139"
 DocStringExtensions = "0.8, 0.9"
 ExponentialUtilities = "1.22"
 FastBroadcast = "0.1.9, 0.2"

src/initdt.jl

@@ -43,12 +43,8 @@
         f(f₀, u0, p, t)
     else
         # TODO: use more caches
-        f₀ = first(get_tmp_cache(integrator))
-        if u0 isa Array && eltype(u0) isa Number
-            fill!(f₀, zero(T))
-        else
-            f₀ = zero.(f₀)
-        end
+        f₀ = DiffEqBase._rate_prototype(u0, t, one(t))
+        f₀ = zero.(f₀)
         f(f₀, u0, p, t)
     end
 

src/solve.jl

@@ -197,20 +197,12 @@ function DiffEqBase.__init(prob::Union{DiffEqBase.AbstractODEProblem,
     dtmax > zero(dtmax) && tdir < 0 && (dtmax *= tdir) # Allow positive dtmax, but auto-convert
     # dtmin is all abs => does not care about sign already.
 
-    if !isdae && isinplace(prob) && u isa AbstractArray && eltype(u) <: Number &&
-       uBottomEltypeNoUnits == uBottomEltype && tType == tTypeNoUnits # Could this be more efficient for other arrays?
-        rate_prototype = recursivecopy(u)
-    elseif prob isa DAEProblem
+    if prob isa DAEProblem
         rate_prototype = prob.du0
     else
-        if (uBottomEltypeNoUnits == uBottomEltype && tType == tTypeNoUnits) ||
-           eltype(u) <: Enum
-            rate_prototype = u
-        else # has units!
-            rate_prototype = u / oneunit(tType)
-        end
+        rate_prototype = DiffEqBase._rate_prototype(u, t, one(t))
     end
-    rateType = typeof(rate_prototype) ## Can be different if united
+    rateType = typeof(rate_prototype) ## Can be different if unitful
 
     if isdae
         if uBottomEltype == uBottomEltypeNoUnits