diff --git a/Project.toml b/Project.toml
index 4386c053c5..25993d7350 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>", "Yingbo Ma <maying
 version = "6.51.2"
 
 [deps]
+ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -32,6 +33,7 @@ RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLNLSolve = "e9a6253c-8580-4d32-9898-8661bb511710"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -45,7 +47,7 @@ TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"
 Adapt = "1.1, 2.0, 3.0"
 ArrayInterface = "6, 7"
 DataStructures = "0.18"
-DiffEqBase = "6.122.0"
+DiffEqBase = "6.125.0"
 DocStringExtensions = "0.8, 0.9"
 ExponentialUtilities = "1.22"
 FastBroadcast = "0.1.9, 0.2"
@@ -67,11 +69,12 @@ Preferences = "1.3"
 RecursiveArrayTools = "2.36"
 Reexport = "0.2, 1.0"
 SciMLBase = "1.90"
+SciMLOperators = "0.2.8"
 SciMLNLSolve = "0.1"
 SimpleNonlinearSolve = "0.1.4"
 SimpleUnPack = "1"
+SparseDiffTools = "2.3"
 PrecompileTools = "1"
-SparseDiffTools = "1.26.2"
 StaticArrayInterface = "1.2"
 StaticArrays = "0.11, 0.12, 1.0"
 TruncatedStacktraces = "1.2"
diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
index 0c89f1cfe6..fa2ba56f97 100644
--- a/src/OrdinaryDiffEq.jl
+++ b/src/OrdinaryDiffEq.jl
@@ -33,7 +33,10 @@ import DiffEqBase: solve!, step!, initialize!, isadaptive
 import DiffEqBase: ODE_DEFAULT_NORM, ODE_DEFAULT_ISOUTOFDOMAIN, ODE_DEFAULT_PROG_MESSAGE,
                    ODE_DEFAULT_UNSTABLE_CHECK
 
-using DiffEqBase: DiffEqArrayOperator, DEFAULT_UPDATE_FUNC
+import SciMLOperators: SciMLOperators, AbstractSciMLOperator, AbstractSciMLScalarOperator,
+                       MatrixOperator, FunctionOperator,
+                       update_coefficients, update_coefficients!, DEFAULT_UPDATE_FUNC,
+                       isconstant
 
 using DiffEqBase: TimeGradientWrapper, UJacobianWrapper, TimeDerivativeWrapper,
                   UDerivativeWrapper
@@ -76,7 +79,7 @@ using FastBroadcast: @.., True, False
 
 using IfElse
 
-using SciMLBase: NoInit, _unwrap_val, AbstractSciMLOperator
+using SciMLBase: NoInit, _unwrap_val
 
 import DiffEqBase: calculate_residuals, calculate_residuals!, unwrap_cache,
                    @tight_loop_macros,
@@ -88,10 +91,12 @@ else
     struct OrdinaryDiffEqTag end
 end
 
-import SparseDiffTools
-import SparseDiffTools: matrix_colors, forwarddiff_color_jacobian!,
+import SparseDiffTools: SparseDiffTools, matrix_colors, forwarddiff_color_jacobian!,
                         forwarddiff_color_jacobian, ForwardColorJacCache,
-                        default_chunk_size, getsize
+                        default_chunk_size, getsize, JacVec
+
+import ADTypes: AbstractADType, AutoFiniteDiff, AutoForwardDiff, AutoReverseDiff,
+                AutoTracker, AutoZygote, AutoEnzyme
 
 import Polyester
 using MacroTools, Adapt
diff --git a/src/alg_utils.jl b/src/alg_utils.jl
index 92cb2c49c3..3f76c0e505 100644
--- a/src/alg_utils.jl
+++ b/src/alg_utils.jl
@@ -17,7 +17,7 @@ function SciMLBase.forwarddiffs_model(alg::Union{OrdinaryDiffEqAdaptiveImplicitA
                                                  DAEAlgorithm,
                                                  OrdinaryDiffEqImplicitAlgorithm,
                                                  ExponentialAlgorithm})
-    alg_autodiff(alg)
+    alg_autodiff(alg) isa AutoForwardDiff
 end
 SciMLBase.forwarddiffs_model_time(alg::RosenbrockAlgorithm) = true
 
@@ -250,7 +250,7 @@ function DiffEqBase.prepare_alg(alg::Union{
                  !(typeof(prob.f.jac_prototype) <: AbstractSciMLOperator)))
             linsolve = LinearSolve.defaultalg(prob.f.jac_prototype, u0)
         else
-            # If mm is a sparse matrix and A is a DiffEqArrayOperator, then let linear
+            # If mm is a sparse matrix and A is a MatrixOperator, then let linear
             # solver choose things later
             linsolve = nothing
         end
@@ -316,20 +316,35 @@ function DiffEqBase.prepare_alg(alg::CompositeAlgorithm, u0, p, prob)
     CompositeAlgorithm(algs, alg.choice_function)
 end
 
-function alg_autodiff(alg::OrdinaryDiffEqAlgorithm)
+# Extract AD type parameter from algorithm, returning as Val to ensure type stability for boolean options.
+function _alg_autodiff(alg::OrdinaryDiffEqAlgorithm)
     error("This algorithm does not have an autodifferentiation option defined.")
 end
-alg_autodiff(alg::OrdinaryDiffEqAdaptiveImplicitAlgorithm{CS, AD}) where {CS, AD} = AD
-alg_autodiff(alg::DAEAlgorithm{CS, AD}) where {CS, AD} = AD
-alg_autodiff(alg::OrdinaryDiffEqImplicitAlgorithm{CS, AD}) where {CS, AD} = AD
-function alg_autodiff(alg::Union{OrdinaryDiffEqExponentialAlgorithm{CS, AD},
-                                 OrdinaryDiffEqAdaptiveExponentialAlgorithm{CS, AD}}) where {
-                                                                                             CS,
-                                                                                             AD
-                                                                                             }
-    AD
+_alg_autodiff(::OrdinaryDiffEqAdaptiveImplicitAlgorithm{CS, AD}) where {CS, AD} = Val{AD}()
+_alg_autodiff(::DAEAlgorithm{CS, AD}) where {CS, AD} = Val{AD}()
+_alg_autodiff(::OrdinaryDiffEqImplicitAlgorithm{CS, AD}) where {CS, AD} = Val{AD}()
+function _alg_autodiff(::Union{OrdinaryDiffEqExponentialAlgorithm{CS, AD},
+                               OrdinaryDiffEqAdaptiveExponentialAlgorithm{CS, AD}
+                              }
+                      ) where {
+                               CS, AD,
+                              }
+    Val{AD}()
 end
 
+function alg_autodiff(alg)
+    autodiff = _alg_autodiff(alg)
+    if autodiff == Val(false)
+        return AutoFiniteDiff()
+    elseif autodiff == Val(true)
+        return AutoForwardDiff()
+    else
+        return _unwrap_val(autodiff)
+    end
+end
+
+# end
+
 # alg_autodiff(alg::CompositeAlgorithm) = alg_autodiff(alg.algs[alg.current_alg])
 get_current_alg_autodiff(alg, cache) = alg_autodiff(alg)
 function get_current_alg_autodiff(alg::CompositeAlgorithm, cache)
diff --git a/src/caches/linear_nonlinear_caches.jl b/src/caches/linear_nonlinear_caches.jl
index ccafa6ac45..92bd228576 100644
--- a/src/caches/linear_nonlinear_caches.jl
+++ b/src/caches/linear_nonlinear_caches.jl
@@ -839,6 +839,19 @@ mutable struct ETD2Fsal{rateType}
     lin::rateType
     nl::rateType
     nlprev::rateType
+
+    function ETD2Fsal(lin, nl, nlprev)
+        if size(lin) == ()
+            # convert to same type if Number or AbstractSciMLScalarOperator
+            T = promote_type(eltype.((lin, nl, nlprev))...)
+
+            lin = convert(T, lin)
+            nl  = convert(T, nl)
+            nlprev  = convert(T, nlprev)
+        end
+
+        new{typeof(lin)}(lin, nl, nlprev)
+    end
 end
 function ETD2Fsal(rate_prototype)
     ETD2Fsal(zero(rate_prototype), zero(rate_prototype), zero(rate_prototype))
diff --git a/src/caches/rosenbrock_caches.jl b/src/caches/rosenbrock_caches.jl
index 1cc59af055..ba1f6da0d8 100644
--- a/src/caches/rosenbrock_caches.jl
+++ b/src/caches/rosenbrock_caches.jl
@@ -152,7 +152,7 @@ function alg_cache(alg::Rosenbrock32, u, rate_prototype, ::Type{uEltypeNoUnits},
                       grad_config, reltol, alg, algebraic_vars)
 end
 
-struct Rosenbrock23ConstantCache{T, TF, UF, JType, WType, F} <: OrdinaryDiffEqConstantCache
+struct Rosenbrock23ConstantCache{T, TF, UF, JType, WType, F, AD} <: OrdinaryDiffEqConstantCache
     c₃₂::T
     d::T
     tf::TF
@@ -160,7 +160,7 @@ struct Rosenbrock23ConstantCache{T, TF, UF, JType, WType, F} <: OrdinaryDiffEqCo
     J::JType
     W::WType
     linsolve::F
-    autodiff::Bool
+    autodiff::AD
 end
 
 function Rosenbrock23ConstantCache(::Type{T}, tf, uf, J, W, linsolve, autodiff) where {T}
@@ -181,7 +181,7 @@ function alg_cache(alg::Rosenbrock23, u, rate_prototype, ::Type{uEltypeNoUnits},
                               alg_autodiff(alg))
 end
 
-struct Rosenbrock32ConstantCache{T, TF, UF, JType, WType, F} <: OrdinaryDiffEqConstantCache
+struct Rosenbrock32ConstantCache{T, TF, UF, JType, WType, F, AD} <: OrdinaryDiffEqConstantCache
     c₃₂::T
     d::T
     tf::TF
@@ -189,7 +189,7 @@ struct Rosenbrock32ConstantCache{T, TF, UF, JType, WType, F} <: OrdinaryDiffEqCo
     J::JType
     W::WType
     linsolve::F
-    autodiff::Bool
+    autodiff::AD
 end
 
 function Rosenbrock32ConstantCache(::Type{T}, tf, uf, J, W, linsolve, autodiff) where {T}
@@ -416,14 +416,14 @@ jac_cache(c::Rosenbrock4Cache) = (c.J, c.W)
 
 ### Rodas methods
 
-struct Rodas4ConstantCache{TF, UF, Tab, JType, WType, F} <: OrdinaryDiffEqConstantCache
+struct Rodas4ConstantCache{TF, UF, Tab, JType, WType, F, AD} <: OrdinaryDiffEqConstantCache
     tf::TF
     uf::UF
     tab::Tab
     J::JType
     W::WType
     linsolve::F
-    autodiff::Bool
+    autodiff::AD
 end
 
 @cache mutable struct Rodas4Cache{uType, rateType, uNoUnitsType, JType, WType, TabType,
diff --git a/src/dense/stiff_addsteps.jl b/src/dense/stiff_addsteps.jl
index e7c497baa8..5e3f68a348 100644
--- a/src/dense/stiff_addsteps.jl
+++ b/src/dense/stiff_addsteps.jl
@@ -7,7 +7,7 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p,
         @unpack tf, uf, d = cache
         γ = dt * d
         tf.u = uprev
-        if cache.autodiff
+        if cache.autodiff isa AutoForwardDiff
             dT = ForwardDiff.derivative(tf, t)
         else
             dT = FiniteDiff.finite_difference_derivative(tf, t, dir = sign(dt))
@@ -180,7 +180,7 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rodas4ConstantCache,
 
         # Time derivative
         tf.u = uprev
-        if cache.autodiff
+        if cache.autodiff isa AutoForwardDiff
             dT = ForwardDiff.derivative(tf, t)
         else
             dT = FiniteDiff.finite_difference_derivative(tf, t, dir = sign(dt))
@@ -566,7 +566,7 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p, cache::Rosenbrock5ConstantCach
 
         # Time derivative
         tf.u = uprev
-        #    if cache.autodiff
+        #    if cache.autodiff isa AutoForwardDiff
         #      dT = ForwardDiff.derivative(tf, t)
         #    else
         dT = FiniteDiff.finite_difference_derivative(tf, t, dir = sign(dt))
diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl
index fd76905816..ff45fca8de 100644
--- a/src/derivative_utils.jl
+++ b/src/derivative_utils.jl
@@ -1,6 +1,6 @@
 const ROSENBROCK_INV_CUTOFF = 7 # https://github.com/SciML/OrdinaryDiffEq.jl/pull/1539
 
-struct StaticWOperator{isinv, T}
+struct StaticWOperator{isinv, T} <: AbstractSciMLOperator{T}
     W::T
     function StaticWOperator(W::T, callinv = true) where {T}
         isinv = size(W, 1) <= ROSENBROCK_INV_CUTOFF
@@ -202,7 +202,7 @@ mutable struct WOperator{IIP, T,
                             transform = false) where {IIP}
         # TODO: there is definitely a missing interface.
         # Tentative interface: `has_concrete` and `concertize(A)`
-        if J isa Union{Number, DiffEqScalar}
+        if J isa Union{Number, ScalarOperator}
             if transform
                 _concrete_form = -mass_matrix / gamma + convert(Number, J)
             else
@@ -210,9 +210,9 @@ mutable struct WOperator{IIP, T,
             end
             _func_cache = nothing
         else
-            AJ = J isa DiffEqArrayOperator ? convert(AbstractMatrix, J) : J
+            AJ = J isa MatrixOperator ? convert(AbstractMatrix, J) : J
             if AJ isa AbstractMatrix
-                mm = mass_matrix isa DiffEqArrayOperator ?
+                mm = mass_matrix isa MatrixOperator ?
                      convert(AbstractMatrix, mass_matrix) : mass_matrix
                 if AJ isa AbstractSparseMatrix
 
@@ -223,7 +223,7 @@ mutable struct WOperator{IIP, T,
                     #
                     # Constant operators never refactorize so always use the correct values there
                     # as well
-                    if gamma == 0 && !(J isa DiffEqArrayOperator && SciMLBase.isconstant(J))
+                    if gamma == 0 && !(J isa MatrixOperator && isconstant(J))
                         # Workaround https://github.com/JuliaSparse/SparseArrays.jl/issues/190
                         # Hopefully `rand()` does not match any value in the array (prob ~ 0, with a check)
                         # Then `one` is required since gamma is zero
@@ -285,7 +285,7 @@ function WOperator{IIP}(f, u, gamma; transform = false) where {IIP}
     J = deepcopy(f.jac_prototype)
     if J isa AbstractMatrix
         @assert DiffEqBase.has_jac(f) "f needs to have an associated jacobian"
-        J = DiffEqArrayOperator(J; update_func = f.jac)
+        J = MatrixOperator(J; update_func! = f.jac)
     end
     return WOperator{IIP}(mass_matrix, gamma, J, u; transform = transform)
 end
@@ -294,17 +294,16 @@ SciMLBase.isinplace(::WOperator{IIP}, i) where {IIP} = IIP
 Base.eltype(W::WOperator) = eltype(W.J)
 
 set_gamma!(W::WOperator, gamma) = (W.gamma = gamma; W)
-function DiffEqBase.update_coefficients!(W::WOperator, u, p, t)
+function SciMLOperators.update_coefficients!(W::WOperator, u, p, t)
     update_coefficients!(W.J, u, p, t)
     update_coefficients!(W.mass_matrix, u, p, t)
-    W.jacvec !== nothing && update_coefficients!(W.jacvec, u, p, t)
+    !isnothing(W.jacvec) && update_coefficients!(W.jacvec, u, p, t)
     W
 end
 
-function DiffEqBase.update_coefficients!(J::SparseDiffTools.JacVec, u, p, t)
-    copyto!(J.x, u)
-    J.f.t = t
-    J.f.p = p
+function SciMLOperators.update_coefficients!(J::UJacobianWrapper, u, p, t)
+    J.p = p
+    J.t = t
 end
 
 function Base.convert(::Type{AbstractMatrix}, W::WOperator{IIP}) where {IIP}
@@ -451,8 +450,7 @@ function do_newJW(integrator, alg, nlsolver, repeat_step)::NTuple{2, Bool}
     isfreshJ = isJcurrent(nlsolver, integrator) && !integrator.u_modified
     iszero(nlsolver.fast_convergence_cutoff) && return isfs && !isfreshJ, isfs
     mm = integrator.f.mass_matrix
-    is_varying_mm = mm isa DiffEqArrayOperator &&
-                    mm.update_func !== SciMLBase.DEFAULT_UPDATE_FUNC
+    is_varying_mm = !isconstant(mm)
     if isfreshJ
         jbad = false
         smallstepchange = true
@@ -675,11 +673,10 @@ function calc_W!(W, integrator, nlsolver::Union{Nothing, AbstractNLSolver}, cach
 
     # calculate W
     if W isa WOperator
-        isnewton(nlsolver) || DiffEqBase.update_coefficients!(W, uprev, p, t) # we will call `update_coefficients!` in NLNewton
+        isnewton(nlsolver) || update_coefficients!(W, uprev, p, t) # we will call `update_coefficients!` in NLNewton
         W.transform = W_transform
         set_gamma!(W, dtgamma)
-        if W.J !== nothing && !(W.J isa SparseDiffTools.JacVec) &&
-           !(W.J isa AbstractSciMLOperator)
+        if W.J !== nothing && !(W.J isa AbstractSciMLOperator)
             islin, isode = islinearfunction(integrator)
             islin ? (J = isode ? f.f : f.f1.f) :
             (new_jac && (calc_J!(W.J, integrator, lcache, next_step)))
@@ -740,7 +737,7 @@ end
         else
             if !isa(J, AbstractSciMLOperator) && (!isnewton(nlsolver) ||
                 nlsolver.cache.W.J isa AbstractSciMLOperator)
-                J = DiffEqArrayOperator(J)
+                J = MatrixOperator(J)
             end
             W = WOperator{false}(mass_matrix, dtgamma, J, uprev, cache.W.jacvec;
                                  transform = W_transform)
@@ -767,7 +764,7 @@ end
         end
     end
     (W isa WOperator && unwrap_alg(integrator, true) isa NewtonAlgorithm) &&
-        (W = DiffEqBase.update_coefficients!(W, uprev, p, t)) # we will call `update_coefficients!` in NLNewton
+        (W = update_coefficients!(W, uprev, p, t)) # we will call `update_coefficients!` in NLNewton
     is_compos && (integrator.eigen_est = isarray ? constvalue(opnorm(J, Inf)) :
                             integrator.opts.internalnorm(J, t))
     return W
@@ -853,8 +850,8 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits},
         # be overridden with concrete_jac.
 
         _f = islin ? (isode ? f.f : f.f1.f) : f
-        jacvec = SparseDiffTools.JacVec(UJacobianWrapper(_f, t, p), copy(u),
-                                        OrdinaryDiffEqTag(), autodiff = alg_autodiff(alg))
+        jacvec = JacVec(UJacobianWrapper(_f, t, p), copy(u), p, t;
+                        autodiff = alg_autodiff(alg), tag = OrdinaryDiffEqTag())
         J = jacvec
         W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec)
 
@@ -869,14 +866,14 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits},
         else
             deepcopy(f.jac_prototype)
         end
-        jacvec = SparseDiffTools.JacVec(UJacobianWrapper(_f, t, p), copy(u),
-                                        OrdinaryDiffEqTag(), autodiff = alg_autodiff(alg))
+        jacvec = JacVec(UJacobianWrapper(_f, t, p), copy(u), p, t;
+                        autodiff = alg_autodiff(alg), tag = OrdinaryDiffEqTag())
         W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec)
 
     elseif islin || (!IIP && DiffEqBase.has_jac(f))
         J = islin ? (isode ? f.f : f.f1.f) : f.jac(uprev, p, t) # unwrap the Jacobian accordingly
         if !isa(J, AbstractSciMLOperator)
-            J = DiffEqArrayOperator(J)
+            J = MatrixOperator(J)
         end
         W = WOperator{IIP}(f.mass_matrix, dt, J, u)
     else
diff --git a/src/derivative_wrappers.jl b/src/derivative_wrappers.jl
index 94317d4222..e63625739e 100644
--- a/src/derivative_wrappers.jl
+++ b/src/derivative_wrappers.jl
@@ -79,8 +79,8 @@ function derivative!(df::AbstractArray{<:Number}, f,
                      x::Union{Number, AbstractArray{<:Number}}, fx::AbstractArray{<:Number},
                      integrator, grad_config)
     alg = unwrap_alg(integrator, true)
-    tmp = length(x) # We calculate derivative for all elements in gradient
-    if alg_autodiff(alg)
+    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)))
         else
@@ -102,7 +102,7 @@ function derivative!(df::AbstractArray{<:Number}, f,
 
         df .= first.(ForwardDiff.partials.(grad_config))
         integrator.stats.nf += 1
-    else
+    elseif alg_autodiff(alg) isa AutoFiniteDiff
         FiniteDiff.finite_difference_gradient!(df, f, x, grad_config,
                                                dir = diffdir(integrator))
         fdtype = alg_difftype(alg)
@@ -113,6 +113,8 @@ function derivative!(df::AbstractArray{<:Number}, f,
             end
         end
         integrator.stats.nf += tmp
+    else
+        error("$alg_autodiff not yet supported in derivative! function")
     end
     nothing
 end
@@ -122,7 +124,7 @@ function derivative(f, x::Union{Number, AbstractArray{<:Number}},
     local d
     tmp = length(x) # We calculate derivative for all elements in gradient
     alg = unwrap_alg(integrator, true)
-    if alg_autodiff(alg)
+    if alg_autodiff(alg) isa AutoForwardDiff
         integrator.stats.nf += 1
         if integrator.iter == 1
             try
@@ -133,7 +135,7 @@ function derivative(f, x::Union{Number, AbstractArray{<:Number}},
         else
             d = ForwardDiff.derivative(f, x)
         end
-    else
+    elseif alg_autodiff(alg) isa AutoFiniteDiff
         d = FiniteDiff.finite_difference_derivative(f, x, alg_difftype(alg),
                                                     dir = diffdir(integrator))
         if alg_difftype(alg) === Val{:central} || alg_difftype(alg) === Val{:forward}
@@ -141,6 +143,8 @@ function derivative(f, x::Union{Number, AbstractArray{<:Number}},
         end
         integrator.stats.nf += tmp
         d
+    else
+        error("$alg_autodiff not yet supported in derivative function")
     end
 end
 
@@ -186,7 +190,7 @@ end
 function jacobian(f, x, integrator)
     alg = unwrap_alg(integrator, true)
     local tmp
-    if alg_autodiff(alg)
+    if alg_autodiff(alg) isa AutoForwardDiff
         if integrator.iter == 1
             try
                 J, tmp = jacobian_autodiff(f, x, integrator.f, alg)
@@ -196,12 +200,14 @@ function jacobian(f, x, integrator)
         else
             J, tmp = jacobian_autodiff(f, x, integrator.f, alg)
         end
-    else
+    elseif alg_autodiff(alg) isa AutoFiniteDiff
         jac_prototype = integrator.f.jac_prototype
         sparsity, colorvec = sparsity_colorvec(integrator.f, x)
         dir = diffdir(integrator)
         J, tmp = jacobian_finitediff(f, x, alg_difftype(alg), dir, colorvec, sparsity,
                                      jac_prototype)
+    else
+        bleh
     end
     integrator.stats.nf += tmp
     J
@@ -222,7 +228,7 @@ function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number},
                    fx::AbstractArray{<:Number}, integrator::DiffEqBase.DEIntegrator,
                    jac_config)
     alg = unwrap_alg(integrator, true)
-    if alg_autodiff(alg)
+    if alg_autodiff(alg) isa AutoForwardDiff 
         if integrator.iter == 1
             try
                 forwarddiff_color_jacobian!(J, f, x, jac_config)
@@ -233,7 +239,7 @@ function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number},
             forwarddiff_color_jacobian!(J, f, x, jac_config)
         end
         integrator.stats.nf += 1
-    else
+    elseif alg_autodiff(alg) isa AutoFiniteDiff 
         isforward = alg_difftype(alg) === Val{:forward}
         if isforward
             forwardcache = get_tmp_cache(integrator, alg, unwrap_cache(integrator, true))[2]
@@ -245,6 +251,8 @@ function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number},
             tmp = jacobian_finitediff!(J, f, x, jac_config, integrator)
         end
         integrator.stats.nf += tmp
+    else 
+        error("$alg_autodiff not yet supported in jacobian! function")
     end
     nothing
 end
@@ -272,7 +280,8 @@ function build_jac_config(alg, f::F1, uf::F2, du1, uprev, u, tmp, du2,
         end
 
         sparsity, colorvec = sparsity_colorvec(f, u)
-        if alg_autodiff(alg)
+
+        if alg_autodiff(alg) isa AutoForwardDiff
             _chunksize = get_chunksize(alg) === Val(0) ? nothing : get_chunksize(alg) # SparseDiffEq uses different convection...
 
             T = if standardtag(alg)
@@ -282,7 +291,7 @@ function build_jac_config(alg, f::F1, uf::F2, du1, uprev, u, tmp, du2,
             end
             jac_config = ForwardColorJacCache(uf, uprev, _chunksize; colorvec = colorvec,
                                               sparsity = sparsity, tag = T)
-        else
+        elseif alg_autodiff(alg) isa AutoFiniteDiff
             if alg_difftype(alg) !== Val{:complex}
                 jac_config = FiniteDiff.JacobianCache(tmp, du1, du2, alg_difftype(alg),
                                                       colorvec = colorvec,
@@ -294,6 +303,8 @@ function build_jac_config(alg, f::F1, uf::F2, du1, uprev, u, tmp, du2,
                                                       colorvec = colorvec,
                                                       sparsity = sparsity)
             end
+        else 
+            error("$alg_autodiff not yet supported in build_jac_config function")
         end
     else
         jac_config = nothing
@@ -343,7 +354,7 @@ end
 
 function build_grad_config(alg, f::F1, tf::F2, du1, t) where {F1, F2}
     if !DiffEqBase.has_tgrad(f)
-        if alg_autodiff(alg)
+        if alg_autodiff(alg) isa AutoForwardDiff
             T = if standardtag(alg)
                 typeof(ForwardDiff.Tag(OrdinaryDiffEqTag(), eltype(du1)))
             else
@@ -362,8 +373,10 @@ function build_grad_config(alg, f::F1, tf::F2, du1, t) where {F1, F2}
                                                                  (ForwardDiff.Partials((one(eltype(du1)),)),)) .*
                                                          false)
             end
-        else
+        elseif alg_autodiff(alg) isa AutoFiniteDiff
             grad_config = FiniteDiff.GradientCache(du1, t, alg_difftype(alg))
+        else 
+            error("$alg_autodiff not yet supported in build_grad_config function")
         end
     else
         grad_config = nothing
diff --git a/src/initialize_dae.jl b/src/initialize_dae.jl
index 95dea03533..23c25105cc 100644
--- a/src/initialize_dae.jl
+++ b/src/initialize_dae.jl
@@ -142,8 +142,8 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation
         failed = nlsolvefail(nlsolver)
         @.. broadcast=false integrator.u=integrator.uprev + z
     else
-        isad = alg_autodiff(integrator.alg)
-        if isad
+        isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff
+        if isAD
             chunk = ForwardDiff.pickchunksize(length(tmp))
             _tmp = PreallocationTools.dualcache(tmp, chunk)
         else
@@ -153,7 +153,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation
         nlequation! = @closure (out, u, p) -> begin
             update_coefficients!(M, u, p, t)
             #M * (u-u0)/dt - f(u,p,t)
-            tmp = isad ? PreallocationTools.get_tmp(_tmp, u) : _tmp
+            tmp = isAD ? PreallocationTools.get_tmp(_tmp, u) : _tmp
             @. tmp = (u - u0) / dt
             mul!(_vec(out), M, _vec(tmp))
             f(tmp, u, p, t)
@@ -161,7 +161,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::ShampineCollocation
             nothing
         end
 
-        nlsolve = default_nlsolve(alg.nlsolve, isinplace, u0, isad)
+        nlsolve = default_nlsolve(alg.nlsolve, isinplace, u0, isAD)
 
         nlfunc = NonlinearFunction(nlequation!;
                                    jac_prototype = f.jac_prototype)
@@ -265,8 +265,8 @@ function _initialize_dae!(integrator, prob::DAEProblem,
     f(resid, integrator.du, u0, p, t)
     integrator.opts.internalnorm(resid, t) <= integrator.opts.abstol && return
 
-    isad = alg_autodiff(integrator.alg)
-    if isad
+    isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff
+    if isAD
         chunk = ForwardDiff.pickchunksize(length(tmp))
         _tmp = PreallocationTools.dualcache(tmp, chunk)
     else
@@ -274,14 +274,14 @@ function _initialize_dae!(integrator, prob::DAEProblem,
     end
 
     nlequation! = @closure (out, u, p) -> begin
-        tmp = isad ? PreallocationTools.get_tmp(_tmp, u) : _tmp
+        tmp = isAD ? PreallocationTools.get_tmp(_tmp, u) : _tmp
         #M * (u-u0)/dt - f(u,p,t)
         @. tmp = (u - u0) / dt
         f(out, tmp, u, p, t)
         nothing
     end
 
-    nlsolve = default_nlsolve(alg.nlsolve, isinplace, u0, isad)
+    nlsolve = default_nlsolve(alg.nlsolve, isinplace, u0, isAD)
 
     nlfunc = NonlinearFunction(nlequation!; jac_prototype = f.jac_prototype)
     nlprob = NonlinearProblem(nlfunc, u0, p)
@@ -372,8 +372,8 @@ function _initialize_dae!(integrator, prob::ODEProblem,
     integrator.opts.internalnorm(tmp, t) <= alg.abstol && return
     alg_u = @view u[algebraic_vars]
 
-    isad = alg_autodiff(integrator.alg)
-    if isad
+    isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff
+    if isAD
         chunk = ForwardDiff.pickchunksize(count(algebraic_vars))
         _tmp = PreallocationTools.dualcache(tmp, chunk)
         _du_tmp = PreallocationTools.dualcache(similar(tmp), chunk)
@@ -382,8 +382,8 @@ function _initialize_dae!(integrator, prob::ODEProblem,
     end
 
     nlequation! = @closure (out, x, p) -> begin
-        uu = isad ? PreallocationTools.get_tmp(_tmp, x) : _tmp
-        du_tmp = isad ? PreallocationTools.get_tmp(_du_tmp, x) : _du_tmp
+        uu = isAD ? PreallocationTools.get_tmp(_tmp, x) : _tmp
+        du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, x) : _du_tmp
         copyto!(uu, integrator.u)
         alg_uu = @view uu[algebraic_vars]
         alg_uu .= x
@@ -394,7 +394,7 @@ function _initialize_dae!(integrator, prob::ODEProblem,
 
     J = algebraic_jacobian(f.jac_prototype, algebraic_eqs, algebraic_vars)
 
-    nlsolve = default_nlsolve(alg.nlsolve, isinplace, u, isad)
+    nlsolve = default_nlsolve(alg.nlsolve, isinplace, u, isAD)
 
     nlfunc = NonlinearFunction(nlequation!; jac_prototype = J)
     nlprob = NonlinearProblem(nlfunc, alg_u, p)
@@ -429,8 +429,8 @@ function _initialize_dae!(integrator, prob::ODEProblem,
 
     integrator.opts.internalnorm(resid, t) <= alg.abstol && return
 
-    isad = alg_autodiff(integrator.alg)
-    if isad
+    isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff
+    if isAD
         chunk = ForwardDiff.pickchunksize(count(algebraic_vars))
         _tmp = PreallocationTools.dualcache(similar(u0), chunk)
     else
@@ -445,7 +445,7 @@ function _initialize_dae!(integrator, prob::ODEProblem,
     end
 
     nlequation = @closure (x, _) -> begin
-        uu = isad ? PreallocationTools.get_tmp(_tmp, x) : _tmp
+        uu = isAD ? PreallocationTools.get_tmp(_tmp, x) : _tmp
         copyto!(uu, integrator.u)
         alg_u = @view uu[algebraic_vars]
         alg_u .= x
@@ -455,7 +455,7 @@ function _initialize_dae!(integrator, prob::ODEProblem,
 
     J = algebraic_jacobian(f.jac_prototype, algebraic_eqs, algebraic_vars)
 
-    nlsolve = default_nlsolve(alg.nlsolve, isinplace, u0, isad)
+    nlsolve = default_nlsolve(alg.nlsolve, isinplace, u0, isAD)
 
     nlfunc = NonlinearFunction(nlequation; jac_prototype = J)
     nlprob = NonlinearProblem(nlfunc, u0[algebraic_vars])
@@ -499,8 +499,8 @@ function _initialize_dae!(integrator, prob::DAEProblem,
         error("differential_vars must be set for DAE initialization to occur. Either set consistent initial conditions, differential_vars, or use a different initialization algorithm.")
     end
 
-    isad = alg_autodiff(integrator.alg)
-    if isad
+    isAD = alg_autodiff(integrator.alg) isa AutoForwardDiff
+    if isAD
         chunk = ForwardDiff.pickchunksize(length(tmp))
         _tmp = PreallocationTools.dualcache(tmp, chunk)
         _du_tmp = PreallocationTools.dualcache(du_tmp, chunk)
@@ -509,8 +509,8 @@ function _initialize_dae!(integrator, prob::DAEProblem,
     end
 
     nlequation! = @closure (out, x, p) -> begin
-        du_tmp = isad ? PreallocationTools.get_tmp(_du_tmp, x) : _du_tmp
-        uu = isad ? PreallocationTools.get_tmp(_tmp, x) : _tmp
+        du_tmp = isAD ? PreallocationTools.get_tmp(_du_tmp, x) : _du_tmp
+        uu = isAD ? PreallocationTools.get_tmp(_tmp, x) : _tmp
 
         @. du_tmp = ifelse(differential_vars, x, du)
         @. uu = ifelse(differential_vars, u, x)
@@ -521,7 +521,7 @@ function _initialize_dae!(integrator, prob::DAEProblem,
     if alg.nlsolve !== nothing
         nlsolve = alg.nlsolve
     else
-        nlsolve = NewtonRaphson(autodiff = isad)
+        nlsolve = NewtonRaphson(autodiff = isAD)
     end
 
     nlfunc = NonlinearFunction(nlequation!; jac_prototype = f.jac_prototype)
diff --git a/src/integrators/integrator_interface.jl b/src/integrators/integrator_interface.jl
index d452119d5d..3cde50a262 100644
--- a/src/integrators/integrator_interface.jl
+++ b/src/integrators/integrator_interface.jl
@@ -234,20 +234,14 @@ function resize_J_W!(cache, integrator, i)
         nf = nlsolve_f(f, integrator.alg)
         islin = f isa Union{ODEFunction, SplitFunction} && islinear(nf.f)
         if !islin
-            if isa(cache.J, AbstractSciMLOperator)
+            if cache.J isa AbstractSciMLOperator
                 resize!(cache.J, i)
             elseif f.jac_prototype !== nothing
                 J = similar(f.jac_prototype, i, i)
-                J = DiffEqArrayOperator(J; update_func = f.jac)
-            elseif cache.J isa SparseDiffTools.JacVec
-                resize!(cache.J.cache1, i)
-                resize!(cache.J.cache2, i)
-                resize!(cache.J.x, i)
+                J = MatrixOperator(J; update_func! = f.jac)
             end
-            if cache.W.jacvec !== nothing
-                resize!(cache.W.jacvec.cache1, i)
-                resize!(cache.W.jacvec.cache2, i)
-                resize!(cache.W.jacvec.x, i)
+            if cache.W.jacvec isa AbstractSciMLOperator
+                resize!(cache.W.jacvec, i)
             end
             cache.W = WOperator{DiffEqBase.isinplace(integrator.sol.prob)}(f.mass_matrix,
                                                                            integrator.dt,
diff --git a/src/misc_utils.jl b/src/misc_utils.jl
index c83d20d1b8..3b7f692f24 100644
--- a/src/misc_utils.jl
+++ b/src/misc_utils.jl
@@ -108,11 +108,13 @@ function dolinsolve(integrator, linsolve; A = nothing, linu = nothing, b = nothi
     # TODO: this ignores the add of the `f` count for add_steps!
     if integrator isa SciMLBase.DEIntegrator && _alg.linsolve !== nothing &&
        !LinearSolve.needs_concrete_A(_alg.linsolve) &&
-       linsolve.A isa WOperator && linsolve.A.J isa SparseDiffTools.JacVec
-        if alg_autodiff(_alg)
+       linsolve.A isa WOperator && linsolve.A.J isa AbstractSciMLOperator 
+        if alg_autodiff(_alg) isa AutoForwardDiff
             integrator.stats.nf += linres.iters
-        else
+        elseif alg_autodiff(_alg) isa AutoFiniteDiff
             integrator.stats.nf += 2 * linres.iters
+        else
+            error("$alg_autodiff not yet supported in dolinsolve function")
         end
     end
 
diff --git a/src/nlsolve/newton.jl b/src/nlsolve/newton.jl
index 49d3f23b52..f00a6f5147 100644
--- a/src/nlsolve/newton.jl
+++ b/src/nlsolve/newton.jl
@@ -169,7 +169,7 @@ end
     b, ustep = _compute_rhs!(nlsolver, integrator, f, z)
 
     # update W
-    if W isa DiffEqBase.AbstractDiffEqLinearOperator
+    if W isa AbstractSciMLOperator
         update_coefficients!(W, ustep, p, tstep)
     end
 
@@ -252,10 +252,10 @@ end
         else
             ustep = @. tmp + γ * z
             if mass_matrix === I
-                ztmp = (dt .* f(ustep, p, tstep) .- z) .* invγdt
+                ztmp = (dt * f(ustep, p, tstep) - z) * invγdt
             else
                 update_coefficients!(mass_matrix, ustep, p, tstep)
-                ztmp = (dt .* f(ustep, p, tstep) .- mass_matrix * z) .* invγdt
+                ztmp = (dt * f(ustep, p, tstep) - mass_matrix * z) * invγdt
             end
         end
     end
diff --git a/src/nlsolve/utils.jl b/src/nlsolve/utils.jl
index f940c261d5..4cb0444dfb 100644
--- a/src/nlsolve/utils.jl
+++ b/src/nlsolve/utils.jl
@@ -61,7 +61,7 @@ function du_alias_or_new(nlsolver::AbstractNLSolver, rate_prototype)
     end
 end
 
-mutable struct DAEResidualJacobianWrapper{AD, F, pType, duType, uType, alphaType, gammaType,
+mutable struct DAEResidualJacobianWrapper{isAD, F, pType, duType, uType, alphaType, gammaType,
                                           tmpType, uprevType, tType} <: Function
     f::F
     p::pType
@@ -73,7 +73,7 @@ mutable struct DAEResidualJacobianWrapper{AD, F, pType, duType, uType, alphaType
     uprev::uprevType
     t::tType
     function DAEResidualJacobianWrapper(alg, f, p, α, invγdt, tmp, uprev, t)
-        isautodiff = alg_autodiff(alg)
+        isautodiff = alg_autodiff(alg) isa AutoForwardDiff
         if isautodiff
             tmp_du = PreallocationTools.dualcache(uprev)
             tmp_u = PreallocationTools.dualcache(uprev)
@@ -87,7 +87,7 @@ mutable struct DAEResidualJacobianWrapper{AD, F, pType, duType, uType, alphaType
     end
 end
 
-is_autodiff(m::DAEResidualJacobianWrapper{AD}) where {AD} = AD
+is_autodiff(m::DAEResidualJacobianWrapper{isAD}) where {isAD} = isAD
 
 function (m::DAEResidualJacobianWrapper)(out, x)
     if is_autodiff(m)
diff --git a/src/perform_step/kencarp_kvaerno_perform_step.jl b/src/perform_step/kencarp_kvaerno_perform_step.jl
index d67824ed60..3087b1b8db 100644
--- a/src/perform_step/kencarp_kvaerno_perform_step.jl
+++ b/src/perform_step/kencarp_kvaerno_perform_step.jl
@@ -206,7 +206,7 @@ end
     if typeof(integrator.f) <: SplitFunction
         # Explicit tableau is not FSAL
         # Make this not compute on repeat
-        z₁ = dt .* f(uprev, p, t)
+        z₁ = dt * f(uprev, p, t)
     else
         # FSAL Step 1
         z₁ = dt * integrator.fsalfirst
diff --git a/src/perform_step/linear_perform_step.jl b/src/perform_step/linear_perform_step.jl
index 0d6546c284..c768bc704f 100644
--- a/src/perform_step/linear_perform_step.jl
+++ b/src/perform_step/linear_perform_step.jl
@@ -795,6 +795,7 @@ function perform_step!(integrator, cache::LinearExponentialCache, repeat_step =
 end
 
 cay!(tmp, A) = mul!(tmp, inv(I - 1 / 2 * A), (I + 1 / 2 * A))
+cay!(tmp, A::AbstractSciMLOperator) = @error "cannot multiply two SciMLOperators with mul!"
 cay(A) = inv(I - 1 / 2 * A) * (I + 1 / 2 * A)
 
 function initialize!(integrator, cache::CayleyEulerConstantCache)
@@ -820,8 +821,10 @@ function perform_step!(integrator, cache::CayleyEulerConstantCache, repeat_step
     end
 
     L = update_coefficients(A, uprev, p, t)
+    L = convert(AbstractMatrix, L)
+
     V = cay(L * dt)
-    u = V * uprev * transpose(V)
+    u = (V * uprev) * transpose(V)
 
     # Update integrator state
     integrator.fsallast = f(u, p, t + dt)
@@ -854,6 +857,7 @@ function perform_step!(integrator, cache::CayleyEulerCache, repeat_step = false)
     end
 
     update_coefficients!(L, uprev, p, t)
+    L = convert(AbstractMatrix, L)
 
     cay!(V, L * dt)
     mul!(tmp, uprev, transpose(V))
diff --git a/test/algconvergence/linear_method_tests.jl b/test/algconvergence/linear_method_tests.jl
index 7e3f44bfa3..7658218e03 100644
--- a/test/algconvergence/linear_method_tests.jl
+++ b/test/algconvergence/linear_method_tests.jl
@@ -2,7 +2,7 @@ using OrdinaryDiffEq, Test, DiffEqDevTools
 using LinearAlgebra, Random
 
 # Linear exponential solvers
-A = DiffEqArrayOperator([2.0 -1.0; -1.0 2.0])
+A = MatrixOperator([2.0 -1.0; -1.0 2.0])
 u0 = ones(2)
 prob = ODEProblem(A, u0, (0.0, 1.0))
 solve(prob, LinearExponential(krylov = :off))
@@ -19,13 +19,13 @@ sol_analytic = exp(1.0 * Matrix(A)) * u0
 @test isapprox(sol4, sol_analytic, rtol = 1e-8)
 
 # u' = A(t)u solvers
-function update_func(A, u, p, t)
+function update_func!(A, u, p, t)
     A[1, 1] = 0
     A[2, 1] = sin(u[1])
     A[1, 2] = -1
     A[2, 2] = 0
 end
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (10, 50.0))
 sol1 = solve(prob, OrdinaryDiffEq.Vern9(), dt = 1 / 4)
 sol2 = solve(prob, OrdinaryDiffEq.RKMK2(), dt = 1 / 4)
@@ -34,7 +34,7 @@ test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
 sim = analyticless_test_convergence(dts, prob, RKMK2(), test_setup)
 @test sim.𝒪est[:l2]≈2 atol=0.2
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (0, 30.0))
 sol1 = solve(prob, OrdinaryDiffEq.Vern9(), dt = 1 / 4)
 sol2 = solve(prob, OrdinaryDiffEq.RKMK4(), dt = 1 / 4)
@@ -43,7 +43,7 @@ test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
 sim = analyticless_test_convergence(dts, prob, RKMK4(), test_setup)
 @test sim.𝒪est[:l2]≈4 atol=0.22
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (0, 30.0))
 sol1 = solve(prob, OrdinaryDiffEq.Vern9(), dt = 1 / 4)
 sol2 = solve(prob, OrdinaryDiffEq.LieRK4(), dt = 1 / 4)
@@ -52,7 +52,7 @@ test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
 sim = analyticless_test_convergence(dts, prob, LieRK4(), test_setup)
 @test sim.𝒪est[:l2]≈5 atol=0.2
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (0, 30.0))
 sol1 = solve(prob, OrdinaryDiffEq.Vern9(), dt = 1 / 4)
 sol2 = solve(prob, OrdinaryDiffEq.CG2(), dt = 1 / 4)
@@ -61,7 +61,7 @@ test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
 sim = analyticless_test_convergence(dts, prob, CG2(), test_setup)
 @test sim.𝒪est[:l2]≈2 atol=0.2
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (0, 20.0))
 sol1 = solve(prob, OrdinaryDiffEq.Vern6(), dt = 1 / 8)
 sol2 = solve(prob, OrdinaryDiffEq.CG3(), dt = 1 / 8)
@@ -70,7 +70,7 @@ test_setup = Dict(:alg => Vern6(), :reltol => 1e-14, :abstol => 1e-14)
 sim = analyticless_test_convergence(dts, prob, CG3(), test_setup)
 @test sim.𝒪est[:l2]≈3 atol=0.2
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (0, 30.0))
 sol1 = solve(prob, Vern9(), dt = 1 / 4)
 sol2 = solve(prob, CG4a(), dt = 1 / 4)
@@ -79,26 +79,26 @@ test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
 sim = analyticless_test_convergence(dts, prob, CG4a(), test_setup)
 @test sim.𝒪est[:l2]≈4 atol=0.2
 
-function update_func(A, u, p, t)
+function update_func!(A, u, p, t)
     A[1, 1] = 0
     A[2, 1] = 1
     A[1, 2] = -2 * (1 - cos(u[2]) - u[2] * sin(u[2]))
     A[2, 2] = 0
 end
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (30, 150.0))
 dts = 1 ./ 2 .^ (7:-1:1)
 test_setup = Dict(:alg => Tsit5(), :reltol => 1e-14, :abstol => 1e-14)
 sim = analyticless_test_convergence(dts, prob, MagnusAdapt4(), test_setup)
 @test sim.𝒪est[:l2]≈4 atol=0.2
 
-function update_func(A, u, p, t)
+function update_func!(A, u, p, t)
     A[1, 1] = cos(t)
     A[2, 1] = sin(t)
     A[1, 2] = -sin(t)
     A[2, 2] = cos(t)
 end
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (0.0, 5.0))
 dts = 1 ./ 2 .^ (10:-1:1)
 sol = solve(prob, OrdinaryDiffEq.MagnusMidpoint(), dt = 1 / 4)
@@ -110,7 +110,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusMidpoint(), test_setup)
 sim = analyticless_test_convergence(dts, prob, MagnusMidpoint(krylov = true), test_setup)
 @test sim.𝒪est[:l2]≈2 atol=0.2
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (0.0, 5.0))
 dts = 1 ./ 2 .^ (10:-1:1)
 sol = solve(prob, OrdinaryDiffEq.MagnusLeapfrog(), dt = 1 / 4)
@@ -122,7 +122,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusLeapfrog(), test_setup)
 sim = analyticless_test_convergence(dts, prob, MagnusLeapfrog(krylov = true), test_setup)
 @test sim.𝒪est[:l2]≈2 atol=0.2
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (0.5, 5.0))
 dts = 1 ./ 2 .^ (10:-1:1)
 sol = solve(prob, OrdinaryDiffEq.LieEuler(), dt = 1 / 4)
@@ -134,7 +134,7 @@ sim = analyticless_test_convergence(dts, prob, LieEuler(), test_setup)
 sim = analyticless_test_convergence(dts, prob, LieEuler(krylov = true), test_setup)
 @test sim.𝒪est[:l2]≈1 atol=0.2
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (1.0, 6.0))
 dts = 1 ./ 2 .^ (10:-1:1)
 sol = solve(prob, MagnusGauss4(), dt = 1 / 4)
@@ -146,7 +146,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusGauss4(), test_setup)
 sim = analyticless_test_convergence(dts, prob, MagnusGauss4(krylov = true), test_setup)
 @test sim.𝒪est[:l2]≈4 atol=0.2
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (1.0, 6.0))
 dts = 1 ./ 2 .^ (4:-1:1)
 test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
@@ -155,7 +155,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusNC6(), test_setup)
 sim = analyticless_test_convergence(dts, prob, MagnusNC6(krylov = true), test_setup)
 @test sim.𝒪est[:l2]≈6 atol=0.2
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (1.0, 6.0))
 sol = solve(prob, MagnusGL6(), dt = 1 / 10)
 dts = 1 ./ 2 .^ (4:-1:1)
@@ -165,7 +165,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusGL6(), test_setup)
 sim = analyticless_test_convergence(dts, prob, MagnusGL6(krylov = true), test_setup)
 @test sim.𝒪est[:l2]≈6 atol=0.3
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (0.0, 100.0))
 dts = 1.775 .^ (5:-1:0)
 test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
@@ -174,7 +174,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusGL8(), test_setup)
 sim = analyticless_test_convergence(dts, prob, MagnusGL8(krylov = true), test_setup)
 @test sim.𝒪est[:l2]≈8 atol=0.2
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (0.0, 100.0))
 dts = 1.773 .^ (5:-1:0)
 test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
@@ -183,7 +183,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusNC8(), test_setup)
 sim = analyticless_test_convergence(dts, prob, MagnusNC8(krylov = true), test_setup)
 @test sim.𝒪est[:l2]≈8 atol=0.2
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (1.0, 6.0))
 dts = 1 ./ 2 .^ (7:-1:1)
 test_setup = Dict(:alg => Vern6(), :reltol => 1e-14, :abstol => 1e-14)
@@ -192,7 +192,7 @@ sim = analyticless_test_convergence(dts, prob, MagnusGL4(), test_setup)
 sim = analyticless_test_convergence(dts, prob, MagnusGL4(krylov = true), test_setup)
 @test sim.𝒪est[:l2]≈4 atol=0.2
 
-A = DiffEqArrayOperator(ones(2, 2), update_func = update_func)
+A = MatrixOperator(ones(2, 2), update_func! = update_func!)
 prob = ODEProblem(A, ones(2), (0, 20.0))
 sol1 = solve(prob, OrdinaryDiffEq.Vern6(), dt = 1 / 8)
 sol2 = solve(prob, OrdinaryDiffEq.CG3(), dt = 1 / 8)
@@ -223,12 +223,16 @@ function B(y::AbstractMatrix)
 end
 
 function update_func(A, u, p, t)
+    B(u)
+end
+
+function update_func!(A, u, p, t)
     A .= B(u)
     return nothing
 end
 
 η = diagm([1.0, 2, 3, 4, 5])
-A = DiffEqArrayOperator(Matrix{eltype(η)}(I(size(η, 1))), update_func = update_func)
+A = MatrixOperator(Matrix{eltype(η)}(I(size(η, 1))), update_func = update_func, update_func! = update_func!)
 dts = 1 ./ 2 .^ (10:-1:2)
 tspan = (0.0, 20.0)
 
diff --git a/test/algconvergence/linear_nonlinear_convergence_tests.jl b/test/algconvergence/linear_nonlinear_convergence_tests.jl
index 06d361fe64..b053431e16 100644
--- a/test/algconvergence/linear_nonlinear_convergence_tests.jl
+++ b/test/algconvergence/linear_nonlinear_convergence_tests.jl
@@ -5,7 +5,7 @@ using OrdinaryDiffEq: alg_order
     println("Caching Out-of-place")
     μ = 1.01
     linnonlin_f2 = (u, p, t) -> μ * u
-    linnonlin_f1 = DiffEqScalar(μ)
+    linnonlin_f1 = ScalarOperator(μ)
     linnonlin_fun = SplitFunction(linnonlin_f1, linnonlin_f2;
                                   analytic = (u0, p, t) -> u0 .* exp.(2μ * t))
     prob = SplitODEProblem(linnonlin_fun, 1 / 2, (0.0, 1.0))
@@ -37,7 +37,7 @@ end
     μ = 1.01
     u0 = rand(2)
     A = [2.0 -1.0; -1.0 2.0]
-    linnonlin_f1 = DiffEqArrayOperator(A)
+    linnonlin_f1 = MatrixOperator(A)
     linnonlin_f2 = (du, u, p, t) -> du .= μ .* u
     linnonlin_fun_iip = SplitFunction(linnonlin_f1, linnonlin_f2;
                                       analytic = (u0, p, t) -> exp((A + μ * I) * t) * u0)
@@ -102,8 +102,8 @@ end
     # Setup nonlinear problem
     A = [-2.0 1.0; 1.0 -2.0]
     f = (du, u, p, t) -> (mul!(du, A, u); du .-= u .^ 3)
-    jac_update = (J, u, p, t) -> (copyto!(J, A); J[1, 1] -= 3u[1]^2; J[2, 2] -= 3u[2]^2)
-    jac_prototype = DiffEqArrayOperator(zeros(2, 2); update_func = jac_update)
+    jac_update! = (J, u, p, t) -> (copyto!(J, A); J[1, 1] -= 3u[1]^2; J[2, 2] -= 3u[2]^2)
+    jac_prototype = MatrixOperator(zeros(2, 2); update_func! = jac_update!)
     fun = ODEFunction(f; jac_prototype = jac_prototype)
     Random.seed!(0)
     u0 = rand(2)
@@ -181,8 +181,8 @@ end
     # Setup nonlinear problem
     A = [-2.0 1.0; 1.0 -2.0]
     f = (du, u, p, t) -> (mul!(du, A, u); du .-= u .^ 3)
-    jac_update = (J, u, p, t) -> (copyto!(J, A); J[1, 1] -= 3u[1]^2; J[2, 2] -= 3u[2]^2)
-    jac_prototype = DiffEqArrayOperator(zeros(2, 2); update_func = jac_update)
+    jac_update! = (J, u, p, t) -> (copyto!(J, A); J[1, 1] -= 3u[1]^2; J[2, 2] -= 3u[2]^2)
+    jac_prototype = MatrixOperator(zeros(2, 2); update_func! = jac_update!)
     fun = ODEFunction(f; jac_prototype = jac_prototype)
     Random.seed!(0)
     u0 = rand(2)
diff --git a/test/algconvergence/linear_nonlinear_krylov_tests.jl b/test/algconvergence/linear_nonlinear_krylov_tests.jl
index 144cbef607..0d46c76149 100644
--- a/test/algconvergence/linear_nonlinear_krylov_tests.jl
+++ b/test/algconvergence/linear_nonlinear_krylov_tests.jl
@@ -8,15 +8,15 @@ let N = 20
     _f = (u, p, t) -> A * u - u .^ 3
     _f_ip = (du, u, p, t) -> (mul!(du, A, u); du .-= u .^ 3)
     _jac = (u, p, t) -> A - 3 * diagm(0 => u .^ 2)
-    _jac_ip = (J, u, p, t) -> begin
+    _jac_ip! = (J, u, p, t) -> begin
         copyto!(J, A)
         @inbounds for i in 1:N
             J[i, i] -= 3 * u[i]^2
         end
     end
     # f = ODEFunction(_f; jac=_jac)
-    # f_ip = ODEFunction(_f_ip; jac=_jac_ip, jac_prototype=zeros(N,N))
-    jac_prototype = DiffEqArrayOperator(zeros(N, N); update_func = _jac_ip)
+    # f_ip = ODEFunction(_f_ip; jac=_jac_ip!, jac_prototype=zeros(N,N))
+    jac_prototype = MatrixOperator(zeros(N, N); update_func! = _jac_ip!)
     f = ODEFunction(_f; jac_prototype = jac_prototype)
     f_ip = ODEFunction(_f_ip; jac_prototype = jac_prototype)
     prob = ODEProblem(f, u0, (0.0, 1.0))
diff --git a/test/algconvergence/ode_rosenbrock_tests.jl b/test/algconvergence/ode_rosenbrock_tests.jl
index 9a130c3fb6..25393948bc 100644
--- a/test/algconvergence/ode_rosenbrock_tests.jl
+++ b/test/algconvergence/ode_rosenbrock_tests.jl
@@ -1,10 +1,11 @@
+using OrdinaryDiffEq, DiffEqDevTools, Test, LinearAlgebra, LinearSolve
+import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear,
+                            prob_ode_bigfloatlinear, prob_ode_bigfloat2Dlinear
+import LinearSolve
+
 @testset "Rosenbrock Tests" begin
     ## Breakout these since no other test of their adaptivity
 
-    using OrdinaryDiffEq, DiffEqDevTools, Test, LinearAlgebra, LinearSolve
-    import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear,
-                              prob_ode_bigfloatlinear, prob_ode_bigfloat2Dlinear
-
     dts = (1 / 2) .^ (6:-1:3)
     testTol = 0.2
 
@@ -435,3 +436,17 @@
     prob = ODEProblem((u, p, t) -> 0.9u, 0.1, (0.0, 1.0))
     @test_nowarn solve(prob, Rosenbrock23(autodiff = false))
 end
+
+@testset "Convergence with time-dependent matrix-free Jacobian" begin
+    time_derivative(du, u, p, t) = (du[1] = t * u[1])
+    time_derivative_analytic(u0, p, t) = u0 * exp(t^2 / 2)
+    ff_time_derivative = ODEFunction(time_derivative, analytic=time_derivative_analytic)
+    prob = ODEProblem(ff_time_derivative, [1.0], (0.0, 1.0))
+
+    dts = (1 / 2) .^ (6:-1:3)
+    testTol = 0.2
+    # Check convergence of Rodas3 with time-dependent matrix-free Jacobian.
+    # Primarily to check that the Jacobian is being updated correctly as t changes.
+    sim = test_convergence(dts, prob, Rodas3(linsolve=LinearSolve.KrylovJL()));
+    @test sim.𝒪est[:final]≈3 atol=testTol
+end
\ No newline at end of file
diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml
index 97aead5433..ac5220e2ba 100644
--- a/test/downstream/Project.toml
+++ b/test/downstream/Project.toml
@@ -4,14 +4,17 @@ DDEProblemLibrary = "f42792ee-6ffc-4e2a-ae83-8ee2f22de800"
 DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
 ODEInterface = "54ca160b-1b9f-5127-a996-1867f4bc2a2c"
 ODEInterfaceDiffEq = "09606e27-ecf5-54fc-bb29-004bd9f985bf"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
-CUDA = "3.12"
+CUDA = "4.2"
 DDEProblemLibrary = "0.1"
-DelayDiffEq = "5.39"
+DelayDiffEq = "5.42"
 ODEInterface = "0.5"
-ODEInterfaceDiffEq = "3.11"
-SciMLSensitivity = "7.11"
-Zygote = "0.6"
\ No newline at end of file
+ODEInterfaceDiffEq = "3.13"
+SciMLSensitivity = "7.30"
+StochasticDiffEq = "6.60.1"
+Zygote = "0.6.61"
diff --git a/test/integrators/ode_event_tests.jl b/test/integrators/ode_event_tests.jl
index 891c0510ce..1b5e410347 100644
--- a/test/integrators/ode_event_tests.jl
+++ b/test/integrators/ode_event_tests.jl
@@ -422,7 +422,7 @@ step!(integrator, 1e-5, true)
     cb = SavingCallback(save_func, saved_values, saveat = t_l)
 
     u0 = normalize(rand(ComplexF64, 100))
-    A = DiffEqArrayOperator(-1im * A)
+    A = MatrixOperator(-1im * A)
     prob = ODEProblem(A, u0, (0, 1.0))
     solve(prob, LinearExponential(), dt = t_l[2] - t_l[1], callback = cb)
     @test length(saved_values.saveval) == length(t_l)
diff --git a/test/interface/ad_tests.jl b/test/interface/ad_tests.jl
index f6d9fe05b0..2af42e88f3 100644
--- a/test/interface/ad_tests.jl
+++ b/test/interface/ad_tests.jl
@@ -329,7 +329,7 @@ fn(x, solver) = objfun(x, prob, data, solver, reltol, abstol)
 
 # test AD with LinearExponential()
 function f(x)
-    K = DiffEqArrayOperator(x)
+    K = MatrixOperator(x)
     u0 = eltype(x).([1.0, 0.0])
     prob = ODEProblem(K, u0, (0.0, 10.0))
     sol = solve(prob, LinearExponential(), tstops = [0.0, 10.0])[2, :]
diff --git a/test/interface/dae_initialize_integration.jl b/test/interface/dae_initialize_integration.jl
index 3315fbdfa4..705d523bee 100644
--- a/test/interface/dae_initialize_integration.jl
+++ b/test/interface/dae_initialize_integration.jl
@@ -37,5 +37,5 @@ p0 = [
 prob = ODEProblem(connected, u0, tspan, p0)
 sol = solve(prob, Rodas5(), initializealg = BrownFullBasicInit())
 @test prob.u0 == sol[1]
-sol = solve(prob, Rodas5(), initiailizealg = ShampineCollocationInit())
+sol = solve(prob, Rodas5(), initializealg = ShampineCollocationInit())
 @test prob.u0 == sol[1]
diff --git a/test/interface/linear_nonlinear_tests.jl b/test/interface/linear_nonlinear_tests.jl
index 8b773b4a39..86cc328668 100644
--- a/test/interface/linear_nonlinear_tests.jl
+++ b/test/interface/linear_nonlinear_tests.jl
@@ -36,35 +36,35 @@ end
 sol = @test_nowarn solve(prob, TRBDF2(autodiff = false));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-                         TRBDF2(autodiff = false, linsolve = IterativeSolversJL_GMRES()));
+                         TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES()));
 @test length(sol.t) < 20
 solref = @test_nowarn solve(prob,
-                            TRBDF2(autodiff = false, linsolve = IterativeSolversJL_GMRES(),
+                            TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(),
                                    smooth_est = false));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-                         TRBDF2(autodiff = false, linsolve = IterativeSolversJL_GMRES(),
+                         TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(),
                                 precs = precsl, smooth_est = false, concrete_jac = true));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-                         TRBDF2(autodiff = false, linsolve = IterativeSolversJL_GMRES(),
+                         TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(),
                                 precs = precsr, smooth_est = false, concrete_jac = true));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-                         TRBDF2(autodiff = false, linsolve = IterativeSolversJL_GMRES(),
+                         TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES(),
                                 precs = precslr, smooth_est = false, concrete_jac = true));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-                         QNDF(autodiff = false, linsolve = IterativeSolversJL_GMRES(),
+                         QNDF(autodiff = false, linsolve = KrylovJL_GMRES(),
                               concrete_jac = true));
 @test length(sol.t) < 25
 sol = @test_nowarn solve(prob,
                          Rosenbrock23(autodiff = false,
-                                      linsolve = IterativeSolversJL_GMRES(),
+                                      linsolve = KrylovJL_GMRES(),
                                       precs = precslr, concrete_jac = true));
 @test length(sol.t) < 20
 sol = @test_nowarn solve(prob,
-                         Rodas4(autodiff = false, linsolve = IterativeSolversJL_GMRES(),
+                         Rodas4(autodiff = false, linsolve = KrylovJL_GMRES(),
                                 precs = precslr, concrete_jac = true));
 @test length(sol.t) < 20
 
diff --git a/test/interface/linear_solver_split_ode_test.jl b/test/interface/linear_solver_split_ode_test.jl
index b4b023aa7f..a3f2ddadbd 100644
--- a/test/interface/linear_solver_split_ode_test.jl
+++ b/test/interface/linear_solver_split_ode_test.jl
@@ -12,8 +12,8 @@ tspan = (0.0, 1.0)
 M1 = 2ones(n) |> Diagonal #|> Array
 M2 = 2ones(n) |> Diagonal #|> Array
 
-f1 = M1 |> DiffEqArrayOperator
-f2 = M2 |> DiffEqArrayOperator
+f1 = M1 |> MatrixOperator
+f2 = M2 |> MatrixOperator
 prob = SplitODEProblem(f1, f2, u0, tspan)
 
 for algname in (:SBDF2,
diff --git a/test/interface/linear_solver_test.jl b/test/interface/linear_solver_test.jl
index f76809670a..374845cfee 100644
--- a/test/interface/linear_solver_test.jl
+++ b/test/interface/linear_solver_test.jl
@@ -5,8 +5,8 @@ using LinearAlgebra, Random
 N = 30
 AA = sprand(MersenneTwister(12), N, N, 0.5)
 mm = sprand(MersenneTwister(123), N, N, 0.5)
-A = DiffEqArrayOperator(AA)
-M = DiffEqArrayOperator(mm'mm)
+A = MatrixOperator(AA)
+M = MatrixOperator(mm'mm)
 u0 = ones(N)
 prob = ODEProblem(ODEFunction(A; mass_matrix = M), u0, (0.0, 1.0))
 
diff --git a/test/interface/mass_matrix_tests.jl b/test/interface/mass_matrix_tests.jl
index 3715910487..0d2093cf5c 100644
--- a/test/interface/mass_matrix_tests.jl
+++ b/test/interface/mass_matrix_tests.jl
@@ -35,7 +35,7 @@ function _norm_dsol(alg, prob, prob2, dt = 1 / 10)
     norm(sol .- sol2)
 end
 
-function update_func1(A, u, p, t)
+function update_func1!(A, u, p, t)
     A[1, 1] = cos(t)
     A[2, 1] = sin(t) * u[1]
     A[3, 1] = t^2
@@ -47,7 +47,13 @@ function update_func1(A, u, p, t)
     A[3, 3] = t * cos(t) + 1
 end
 
-function update_func2(A, u, p, t)
+function update_func1(A, u, p, t)
+    A = copy(A)
+    update_func1!(A, u, p, t)
+    A
+end
+
+function update_func2!(A, u, p, t)
     A[1, 1] = cos(t)
     A[2, 1] = sin(t)
     A[3, 1] = t^2
@@ -59,10 +65,16 @@ function update_func2(A, u, p, t)
     A[3, 3] = t * cos(t) + 1
 end
 
+function update_func2(A, u, p, t)
+    A = copy(A)
+    update_func2!(A, u, p, t)
+    A
+end
+
 almost_I = Matrix{Float64}(1.01I, 3, 3)
 mm_A = Float64[-2 1 4; 4 -2 1; 2 1 3]
-dependent_M1 = DiffEqArrayOperator(ones(3, 3), update_func = update_func1)
-dependent_M2 = DiffEqArrayOperator(ones(3, 3), update_func = update_func2)
+dependent_M1 = MatrixOperator(ones(3, 3), update_func = update_func1, update_func! = update_func1!)
+dependent_M2 = MatrixOperator(ones(3, 3), update_func = update_func2, update_func! = update_func2!)
 
 @testset "Mass Matrix Accuracy Tests" for mm in (almost_I, mm_A)
     # test each method for exactness
diff --git a/test/interface/nojac.jl b/test/interface/nojac.jl
index 86800bf824..0696e7101f 100644
--- a/test/interface/nojac.jl
+++ b/test/interface/nojac.jl
@@ -5,6 +5,9 @@ const xyd_brusselator = range(0, stop = 1, length = N)
 brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0
 limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a
 function brusselator_2d_loop(du, u, p, t)
+    u = reshape(u, N, N, 2)
+    du = reshape(du, N, N, 2)
+
     A, B, alpha, dx = p
     alpha = alpha / dx^2
     @inbounds for I in CartesianIndices((N, N))
@@ -20,6 +23,7 @@ function brusselator_2d_loop(du, u, p, t)
                        4u[i, j, 2]) +
                       A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2]
     end
+    vec(du)
 end
 p = (3.4, 1.0, 10.0, step(xyd_brusselator))
 
@@ -34,7 +38,7 @@ function init_brusselator_2d(xyd)
     end
     u
 end
-u0 = init_brusselator_2d(xyd_brusselator)
+u0 = init_brusselator_2d(xyd_brusselator) |> vec
 prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
                                      u0, (0.0, 5.0), p)
 
diff --git a/test/interface/nonfulldiagonal_sparse.jl b/test/interface/nonfulldiagonal_sparse.jl
index 6f15ef050e..27018a334b 100644
--- a/test/interface/nonfulldiagonal_sparse.jl
+++ b/test/interface/nonfulldiagonal_sparse.jl
@@ -19,7 +19,7 @@ function enclosethetimedifferential(parameters::NamedTuple)::Function
         lower[1] = -1.0
         upper[end] = 1.0
         M = hcat(lower, sparse(diagm(-1 => du, 0 => diag, 1 => du2)), upper)
-        DiffEqArrayOperator(1 / dx * M)
+        MatrixOperator(1 / dx * M)
     end
 
     function second_deriv(N)
@@ -32,7 +32,7 @@ function enclosethetimedifferential(parameters::NamedTuple)::Function
         lower[1] = 1.0
         upper[end] = 1.0
         M = hcat(lower, sparse(diagm(-1 => du, 0 => diag, 1 => du2)), upper)
-        DiffEqArrayOperator(1 / dx^2 * M)
+        MatrixOperator(1 / dx^2 * M)
     end
 
     function extender(N)
@@ -45,7 +45,7 @@ function enclosethetimedifferential(parameters::NamedTuple)::Function
         M = vcat(transpose(lower),
                  sparse(diagm(diag)),
                  transpose(upper))
-        DiffEqArrayOperator(1 / dx^2 * M)
+        MatrixOperator(1 / dx^2 * M)
     end
 
     bc_handler = extender(N)
diff --git a/test/interface/norecompile.jl b/test/interface/norecompile.jl
index 1443b76265..48f9ff95cc 100644
--- a/test/interface/norecompile.jl
+++ b/test/interface/norecompile.jl
@@ -14,14 +14,19 @@ function lorenz(du, u, p, t)
 end
 lorenzprob = ODEProblem(lorenz, [1.0; 0.0; 0.0], (0.0, 1.0), Float64[])
 
-t1 = @elapsed sol = solve(lorenzprob, Rosenbrock23())
-t2 = @elapsed sol = solve(lorenzprob, Rosenbrock23(autodiff = false))
+t1 = @elapsed sol1 = solve(lorenzprob, Rosenbrock23())
+t2 = @elapsed sol2 = solve(lorenzprob, Rosenbrock23(autodiff = false))
 
 lorenzprob2 = ODEProblem{true, SciMLBase.FullSpecialize}(lorenz, [1.0; 0.0; 0.0],
                                                          (0.0, 1.0), Float64[])
 
-t3 = @elapsed sol = solve(lorenzprob2, Rosenbrock23())
-t4 = @elapsed sol = solve(lorenzprob2, Rosenbrock23(autodiff = false))
+t3 = @elapsed sol3 = solve(lorenzprob2, Rosenbrock23())
+t4 = @elapsed sol4 = solve(lorenzprob2, Rosenbrock23(autodiff = false))
+
+@test sol1.retcode === ReturnCode.Success
+@test sol2.retcode === ReturnCode.Success
+@test sol3.retcode === ReturnCode.Success
+@test sol4.retcode === ReturnCode.Success
 
 if VERSION >= v"1.8"
     @test t1 < t3
diff --git a/test/interface/preconditioners.jl b/test/interface/preconditioners.jl
index 5d6180c56c..8cabf7de9a 100644
--- a/test/interface/preconditioners.jl
+++ b/test/interface/preconditioners.jl
@@ -13,6 +13,10 @@ limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a
 const iter = Ref(0)
 function brusselator_2d_loop(du, u, p, t)
     global iter[] += 1
+
+    u = reshape(u, N, N, 2)
+    du = reshape(du, N, N, 2)
+
     A, B, alpha, dx = p
     alpha = alpha / dx^2
     @inbounds for I in CartesianIndices((N, N))
@@ -28,6 +32,8 @@ function brusselator_2d_loop(du, u, p, t)
                        4u[i, j, 2]) +
                       A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2]
     end
+
+    vec(du)
 end
 p = (3.4, 1.0, 10.0, step(xyd_brusselator))
 
@@ -42,9 +48,8 @@ function init_brusselator_2d(xyd)
     end
     u
 end
-u0 = init_brusselator_2d(xyd_brusselator)
-prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
-                                     u0, (0.0, 11.5), p)
+u0 = init_brusselator_2d(xyd_brusselator) |> vec
+prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop, u0, (0.0, 11.5), p)
 
 du0 = copy(u0)
 jac = ModelingToolkit.Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p,
@@ -109,6 +114,11 @@ sol4 = solve(prob_ode_brusselator_2d_sparse,
 iter4 = iter[];
 iter[] = 0;
 
+@test sol1.retcode === ReturnCode.Success
+@test sol2.retcode === ReturnCode.Success
+@test sol3.retcode === ReturnCode.Success
+@test sol4.retcode === ReturnCode.Success
+
 @test iter2 < iter1
 @test iter3 < iter1
 @test iter4 < iter1
@@ -133,6 +143,11 @@ sol4 = solve(prob_ode_brusselator_2d_sparse,
 iter4 = iter[];
 iter[] = 0;
 
+@test sol1.retcode === ReturnCode.Success
+@test sol2.retcode === ReturnCode.Success
+@test sol3.retcode === ReturnCode.Success
+@test sol4.retcode === ReturnCode.Success
+
 @test iter2 < iter1
 @test iter3 < iter1
 @test iter4 < iter1
@@ -157,6 +172,11 @@ sol4 = solve(prob_ode_brusselator_2d_sparse,
 iter4 = iter[];
 iter[] = 0;
 
+@test sol1.retcode === ReturnCode.Success
+@test sol2.retcode === ReturnCode.Success
+@test sol3.retcode === ReturnCode.Success
+@test sol4.retcode === ReturnCode.Success
+
 @test iter2 < iter1
 @test iter3 < iter1
 @test iter4 < iter1
@@ -181,6 +201,11 @@ sol4 = solve(prob_ode_brusselator_2d_sparse,
 iter4 = iter[];
 iter[] = 0;
 
+@test sol1.retcode === ReturnCode.Success
+@test sol2.retcode === ReturnCode.Success
+@test sol3.retcode === ReturnCode.Success
+@test sol4.retcode === ReturnCode.Success
+
 @test iter2 < iter1
 @test iter3 < iter1
 @test iter4 < iter1
@@ -205,6 +230,11 @@ sol4 = solve(prob_ode_brusselator_2d_sparse,
 iter4 = iter[];
 iter[] = 0;
 
+@test sol1.retcode === ReturnCode.Success
+@test sol2.retcode === ReturnCode.Success
+@test sol3.retcode === ReturnCode.Success
+@test sol4.retcode === ReturnCode.Success
+
 @test iter2 < iter1
 @test iter3 < iter1
 @test iter4 < iter1
@@ -229,6 +259,11 @@ sol4 = solve(prob_ode_brusselator_2d_sparse,
 iter4 = iter[];
 iter[] = 0;
 
+@test sol1.retcode === ReturnCode.Success
+@test sol2.retcode === ReturnCode.Success
+@test sol3.retcode === ReturnCode.Success
+@test sol4.retcode === ReturnCode.Success
+
 @test iter2 < iter1
 @test iter3 < iter1
 @test iter4 < iter1
diff --git a/test/interface/utility_tests.jl b/test/interface/utility_tests.jl
index 934892203d..44fc551812 100644
--- a/test/interface/utility_tests.jl
+++ b/test/interface/utility_tests.jl
@@ -24,7 +24,7 @@ using OrdinaryDiffEq, LinearAlgebra, SparseArrays, Random, Test, LinearSolve
     # In-place
     fun = ODEFunction((du, u, p, t) -> mul!(du, A, u);
                       mass_matrix = mm,
-                      jac_prototype = DiffEqArrayOperator(A))
+                      jac_prototype = MatrixOperator(A))
     integrator = init(ODEProblem(fun, u0, tspan), ImplicitEuler(); adaptive = false,
                       dt = dt)
     calc_W!(integrator.cache.nlsolver.cache.W, integrator, integrator.cache.nlsolver,
@@ -56,10 +56,7 @@ end
     fun2 = ODEFunction(_f; mass_matrix = mm, jac = (u, p, t) -> t * A)
     fun1_ip = ODEFunction(_f_ip; mass_matrix = mm)
     fun2_ip = ODEFunction(_f_ip; mass_matrix = mm,
-                          jac_prototype = DiffEqArrayOperator(similar(A);
-                                                              update_func = (J, u, p, t) -> (J .= t .*
-                                                                                                  A;
-                                                                                             J)))
+                          jac_prototype = MatrixOperator(similar(A); update_func! = (J, u, p, t) -> J .= t .* A))
 
     for Alg in [ImplicitEuler, Rosenbrock23, Rodas5]
         println(Alg)