fix the solve

src/derivative_utils.jl

@@ -294,7 +294,7 @@ SciMLBase.isinplace(::WOperator{IIP}, i) where {IIP} = IIP
 Base.eltype(W::WOperator) = eltype(W.J)
 
 # In WOperator update_coefficients!, accept both missing u/p/t and missing dtgamma/transform and don't update them in that case.
-# This helps support partial updating logic used with Newton solvers. 
+# This helps support partial updating logic used with Newton solvers.
 function SciMLOperators.update_coefficients!(W::WOperator,
     u = nothing,
     p = nothing,

src/derivative_wrappers.jl

@@ -82,13 +82,12 @@ function derivative!(df::AbstractArray{<:Number}, f,
     tmp = length(x) # We calculate derivtive for all elements in gradient
     if alg_autodiff(alg) isa AutoForwardDiff
         T = if standardtag(alg)
-            typeof(ForwardDiff.Tag(OrdinaryDiffEqTag(), eltype(df)))
+            typeof(ForwardDiff.Tag(OrdinaryDiffEqTag(), eltype(x)))
         else
-            typeof(ForwardDiff.Tag(f, eltype(df)))
+            typeof(ForwardDiff.Tag(f, eltype(x)))
         end
-
-        xdual = Dual{T, eltype(df), 1}(convert(eltype(df), x),
-            ForwardDiff.Partials((one(eltype(df)),)))
+        xdual = Dual{T, eltype(x), 1}(x,
+            ForwardDiff.Partials((one(eltype(x)),)))
 
         if integrator.iter == 1
             try

src/initdt.jl

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