diff --git a/Project.toml b/Project.toml
index a06a02b29e..dcb1503466 100644
--- a/Project.toml
+++ b/Project.toml
@@ -85,7 +85,7 @@ RecipesBase = "1.3.4"
 RecursiveArrayTools = "3.27.2"
 Reexport = "1"
 RuntimeGeneratedFunctions = "0.5.12"
-SciMLOperators = "0.3.7"
+SciMLOperators = "0.3.13"
 SciMLStructures = "1.1"
 StableRNGs = "1.0"
 StaticArrays = "1.7"
diff --git a/src/SciMLBase.jl b/src/SciMLBase.jl
index 526e7589cc..92bc2d8c3a 100644
--- a/src/SciMLBase.jl
+++ b/src/SciMLBase.jl
@@ -575,21 +575,6 @@ $(TYPEDEF)
 abstract type AbstractSensitivitySolution{T, N, S} <: AbstractTimeseriesSolution{T, N, S} end
 
 # Misc
-# TODO - deprecate AbstractDiffEqOperator family
-"""
-$(TYPEDEF)
-"""
-abstract type AbstractDiffEqOperator{T} <: AbstractSciMLOperator{T} end
-
-"""
-$(TYPEDEF)
-"""
-abstract type AbstractDiffEqLinearOperator{T} <: AbstractDiffEqOperator{T} end
-
-"""
-$(TYPEDEF)
-"""
-abstract type AbstractDiffEqCompositeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 
 """
 $(TYPEDEF)
@@ -720,10 +705,6 @@ $(TYPEDEF)
 abstract type AbstractParameterizedFunction{iip} <: AbstractODEFunction{iip} end
 
 include("retcodes.jl")
-include("operators/operators.jl")
-include("operators/basic_operators.jl")
-include("operators/diffeq_operator.jl")
-include("operators/common_defaults.jl")
 include("symbolic_utils.jl")
 include("performance_warnings.jl")
 
@@ -853,10 +834,6 @@ export EnsembleAnalysis, EnsembleSummary
 
 export tuples, intervals, TimeChoiceIterator
 
-export AffineDiffEqOperator, DiffEqScaledOperator
-
-export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity
-
 export step!, deleteat!, addat!, get_tmp_cache,
        full_cache, user_cache, u_cache, du_cache,
        rand_cache, ratenoise_cache,
diff --git a/src/operators/basic_operators.jl b/src/operators/basic_operators.jl
deleted file mode 100644
index 450206e30a..0000000000
--- a/src/operators/basic_operators.jl
+++ /dev/null
@@ -1,206 +0,0 @@
-"""
-$(TYPEDEF)
-"""
-struct DiffEqIdentity{T, N} <: AbstractDiffEqLinearOperator{T} end
-
-DiffEqIdentity(u) = DiffEqIdentity{eltype(u), length(u)}()
-Base.size(::DiffEqIdentity{T, N}) where {T, N} = (N, N)
-Base.size(::DiffEqIdentity{T, N}, m::Integer) where {T, N} = (m == 1 || m == 2) ? N : 1
-LinearAlgebra.opnorm(::DiffEqIdentity{T, N}, p::Real = 2) where {T, N} = one(T)
-function Base.convert(::Type{AbstractMatrix}, ::DiffEqIdentity{T, N}) where {T, N}
-    LinearAlgebra.Diagonal(ones(T, N))
-end
-
-for op in (:*, :/, :\)
-    @eval Base.$op(::DiffEqIdentity{T, N}, x::AbstractVecOrMat) where {T, N} = $op(I, x)
-    @eval Base.$op(::DiffEqIdentity{T, N}, x::AbstractArray) where {T, N} = $op(I, x)
-    @eval Base.$op(x::AbstractVecOrMat, ::DiffEqIdentity{T, N}) where {T, N} = $op(x, I)
-    @eval Base.$op(x::AbstractArray, ::DiffEqIdentity{T, N}) where {T, N} = $op(x, I)
-end
-LinearAlgebra.mul!(Y::AbstractVecOrMat, ::DiffEqIdentity, B::AbstractVecOrMat) = Y .= B
-LinearAlgebra.ldiv!(Y::AbstractVecOrMat, ::DiffEqIdentity, B::AbstractVecOrMat) = Y .= B
-
-LinearAlgebra.mul!(Y::AbstractArray, ::DiffEqIdentity, B::AbstractArray) = Y .= B
-LinearAlgebra.ldiv!(Y::AbstractArray, ::DiffEqIdentity, B::AbstractArray) = Y .= B
-
-for pred in (:isreal, :issymmetric, :ishermitian, :isposdef)
-    @eval LinearAlgebra.$pred(::DiffEqIdentity) = true
-end
-
-"""
-    DiffEqScalar(val[; update_func])
-
-Represents a time-dependent scalar/scaling operator. The update function
-is called by `update_coefficients!` and is assumed to have the following
-signature:
-
-    update_func(oldval,u,p,t) -> newval
-"""
-mutable struct DiffEqScalar{T <: Number, F} <: AbstractDiffEqLinearOperator{T}
-    val::T
-    update_func::F
-    function DiffEqScalar(val::T; update_func = DEFAULT_UPDATE_FUNC) where {T}
-        new{T, typeof(update_func)}(val, update_func)
-    end
-end
-
-Base.convert(::Type{Number}, α::DiffEqScalar) = α.val
-Base.convert(::Type{DiffEqScalar}, α::Number) = DiffEqScalar(α)
-Base.size(::DiffEqScalar) = ()
-Base.size(::DiffEqScalar, ::Integer) = 1
-update_coefficients!(α::DiffEqScalar, u, p, t) = (α.val = α.update_func(α.val, u, p, t);
-α)
-isconstant(α::DiffEqScalar) = α.update_func == DEFAULT_UPDATE_FUNC
-
-for op in (:*, :/, :\)
-    for T in (
-        ### added in https://github.com/SciML/SciMLBase.jl/pull/377
-        :AbstractVecOrMat,
-        ###
-        :AbstractArray,
-        :Number)
-        @eval Base.$op(α::DiffEqScalar, x::$T) = $op(α.val, x)
-        @eval Base.$op(x::$T, α::DiffEqScalar) = $op(x, α.val)
-    end
-    @eval Base.$op(x::DiffEqScalar, y::DiffEqScalar) = $op(x.val, y.val)
-end
-
-for op in (:-, :+)
-    @eval Base.$op(α::DiffEqScalar, x::Number) = $op(α.val, x)
-    @eval Base.$op(x::Number, α::DiffEqScalar) = $op(x, α.val)
-    @eval Base.$op(x::DiffEqScalar, y::DiffEqScalar) = $op(x.val, y.val)
-end
-
-LinearAlgebra.lmul!(α::DiffEqScalar, B::AbstractArray) = lmul!(α.val, B)
-LinearAlgebra.rmul!(B::AbstractArray, α::DiffEqScalar) = rmul!(B, α.val)
-LinearAlgebra.mul!(Y::AbstractArray, α::DiffEqScalar, B::AbstractArray) = mul!(Y, α.val, B)
-function LinearAlgebra.axpy!(α::DiffEqScalar, X::AbstractArray, Y::AbstractArray)
-    axpy!(α.val, X, Y)
-end
-Base.abs(α::DiffEqScalar) = abs(α.val)
-
-"""
-    DiffEqArrayOperator(A[; update_func])
-
-Represents a time-dependent linear operator given by an AbstractMatrix. The
-update function is called by `update_coefficients!` and is assumed to have
-the following signature:
-
-    update_func(A::AbstractMatrix,u,p,t) -> [modifies A]
-"""
-struct DiffEqArrayOperator{T, AType <: AbstractMatrix{T}, F} <:
-       AbstractDiffEqLinearOperator{T}
-    A::AType
-    update_func::F
-    function DiffEqArrayOperator(A::AType; update_func = DEFAULT_UPDATE_FUNC) where {AType}
-        new{eltype(A), AType, typeof(update_func)}(A, update_func)
-    end
-end
-
-has_adjoint(::DiffEqArrayOperator) = true
-update_coefficients!(L::DiffEqArrayOperator, u, p, t) = (L.update_func(L.A, u, p, t); L)
-isconstant(L::DiffEqArrayOperator) = L.update_func == DEFAULT_UPDATE_FUNC
-function Base.similar(L::DiffEqArrayOperator, ::Type{T}, dims::Dims) where {T}
-    similar(L.A, T, dims)
-end
-function Base.adjoint(L::DiffEqArrayOperator)
-    DiffEqArrayOperator(L.A'; update_func = (A, u, p, t) -> L.update_func(L.A, u, p, t)')
-end
-
-# propagate_inbounds here for the getindex fallback
-Base.@propagate_inbounds Base.convert(::Type{AbstractMatrix}, L::DiffEqArrayOperator) = L.A
-Base.@propagate_inbounds Base.setindex!(L::DiffEqArrayOperator, v, i::Int) = (L.A[i] = v)
-Base.@propagate_inbounds function Base.setindex!(L::DiffEqArrayOperator, v,
-        I::Vararg{Int, N}) where {N}
-    (L.A[I...] = v)
-end
-
-Base.eachcol(L::DiffEqArrayOperator) = eachcol(L.A)
-Base.eachrow(L::DiffEqArrayOperator) = eachrow(L.A)
-Base.length(L::DiffEqArrayOperator) = length(L.A)
-Base.iterate(L::DiffEqArrayOperator, args...) = iterate(L.A, args...)
-Base.axes(L::DiffEqArrayOperator) = axes(L.A)
-Base.eachindex(L::DiffEqArrayOperator) = eachindex(L.A)
-function Base.IndexStyle(::Type{<:DiffEqArrayOperator{T, AType}}) where {T, AType}
-    Base.IndexStyle(AType)
-end
-Base.copyto!(L::DiffEqArrayOperator, rhs) = (copyto!(L.A, rhs); L)
-Base.Broadcast.broadcastable(L::DiffEqArrayOperator) = L
-Base.ndims(::Type{<:DiffEqArrayOperator{T, AType}}) where {T, AType} = ndims(AType)
-ArrayInterface.issingular(L::DiffEqArrayOperator) = ArrayInterface.issingular(L.A)
-function Base.copy(L::DiffEqArrayOperator)
-    DiffEqArrayOperator(copy(L.A); update_func = L.update_func)
-end
-
-const AdjointFact = isdefined(LinearAlgebra, :AdjointFactorization) ?
-                    LinearAlgebra.AdjointFactorization : Adjoint
-const TransposeFact = isdefined(LinearAlgebra, :TransposeFactorization) ?
-                      LinearAlgebra.TransposeFactorization : Transpose
-
-"""
-    FactorizedDiffEqArrayOperator(F)
-
-Like DiffEqArrayOperator, but stores a Factorization instead.
-
-Supports left division and `ldiv!` when applied to an array.
-"""
-struct FactorizedDiffEqArrayOperator{T <: Number,
-    FType <: Union{
-        Factorization{T}, Diagonal{T}, Bidiagonal{T},
-        AdjointFact{T, <:Factorization{T}},
-        TransposeFact{T, <:Factorization{T}}
-    }
-} <: AbstractDiffEqLinearOperator{T}
-    F::FType
-end
-
-function Base.convert(::Type{AbstractMatrix},
-        L::FactorizedDiffEqArrayOperator{<:Any,
-            <:Union{Factorization, AbstractMatrix
-            }})
-    convert(AbstractMatrix, L.F)
-end
-function Base.convert(::Type{AbstractMatrix},
-        L::FactorizedDiffEqArrayOperator{<:Any, <:Union{Adjoint, AdjointFact}
-        })
-    adjoint(convert(AbstractMatrix, adjoint(L.F)))
-end
-function Base.convert(::Type{AbstractMatrix},
-        L::FactorizedDiffEqArrayOperator{<:Any,
-            <:Union{Transpose, TransposeFact}})
-    transpose(convert(AbstractMatrix, transpose(L.F)))
-end
-
-function Base.Matrix(L::FactorizedDiffEqArrayOperator{<:Any,
-        <:Union{Factorization, AbstractMatrix
-        }})
-    Matrix(L.F)
-end
-function Base.Matrix(L::FactorizedDiffEqArrayOperator{<:Any, <:Union{Adjoint, AdjointFact}})
-    adjoint(Matrix(adjoint(L.F)))
-end
-function Base.Matrix(L::FactorizedDiffEqArrayOperator{<:Any,
-        <:Union{Transpose, TransposeFact}})
-    transpose(Matrix(transpose(L.F)))
-end
-
-Base.adjoint(L::FactorizedDiffEqArrayOperator) = FactorizedDiffEqArrayOperator(L.F')
-Base.size(L::FactorizedDiffEqArrayOperator, args...) = size(L.F, args...)
-function LinearAlgebra.ldiv!(Y::AbstractVecOrMat, L::FactorizedDiffEqArrayOperator,
-        B::AbstractVecOrMat)
-    ldiv!(Y, L.F, B)
-end
-LinearAlgebra.ldiv!(L::FactorizedDiffEqArrayOperator, B::AbstractVecOrMat) = ldiv!(L.F, B)
-Base.:\(L::FactorizedDiffEqArrayOperator, x::AbstractVecOrMat) = L.F \ x
-LinearAlgebra.issuccess(L::FactorizedDiffEqArrayOperator) = issuccess(L.F)
-
-LinearAlgebra.ldiv!(y, L::FactorizedDiffEqArrayOperator, x) = ldiv!(y, L.F, x)
-#isconstant(::FactorizedDiffEqArrayOperator) = true
-has_ldiv(::FactorizedDiffEqArrayOperator) = true
-has_ldiv!(::FactorizedDiffEqArrayOperator) = true
-
-# The (u,p,t) and (du,u,p,t) interface
-for T in [DiffEqScalar, DiffEqArrayOperator, FactorizedDiffEqArrayOperator, DiffEqIdentity]
-    (L::T)(u, p, t) = (update_coefficients!(L, u, p, t); L * u)
-    (L::T)(du, u, p, t) = (update_coefficients!(L, u, p, t); mul!(du, L, u))
-end
diff --git a/src/operators/common_defaults.jl b/src/operators/common_defaults.jl
deleted file mode 100644
index 632b147828..0000000000
--- a/src/operators/common_defaults.jl
+++ /dev/null
@@ -1,91 +0,0 @@
-update_coefficients!(L::AbstractDiffEqLinearOperator, u, p, t) = L
-
-# Routines that use the AbstractMatrix representation
-function Base.convert(::Type{AbstractArray}, L::AbstractDiffEqLinearOperator)
-    convert(AbstractMatrix, L)
-end
-function Base.size(L::AbstractDiffEqLinearOperator, i::Integer)
-    size(convert(AbstractMatrix, L), i)
-end
-function Base.size(L::AbstractDiffEqLinearOperator, args...)
-    size(convert(AbstractMatrix, L), args...)
-end
-function LinearAlgebra.opnorm(L::AbstractDiffEqLinearOperator, p::Real = 2)
-    opnorm(convert(AbstractMatrix, L), p)
-end
-Base.@propagate_inbounds function Base.getindex(L::AbstractDiffEqLinearOperator,
-        I::Vararg{Any, N}) where {N}
-    convert(AbstractMatrix, L)[I...]
-end
-function Base.getindex(L::AbstractDiffEqLinearOperator, I::Vararg{Int, N}) where {N}
-    convert(AbstractMatrix, L)[I...]
-end
-for op in (:*, :/, :\)
-    ### added in https://github.com/SciML/SciMLBase.jl/pull/377
-    @eval function Base.$op(L::AbstractDiffEqLinearOperator, x::AbstractVecOrMat)
-        $op(convert(AbstractMatrix, L), x)
-    end
-    @eval function Base.$op(x::AbstractVecOrMat, L::AbstractDiffEqLinearOperator)
-        $op(x, convert(AbstractMatrix, L))
-    end
-    ###
-    @eval function Base.$op(L::AbstractDiffEqLinearOperator, x::AbstractArray)
-        $op(convert(AbstractMatrix, L), x)
-    end
-    @eval function Base.$op(x::AbstractArray, L::AbstractDiffEqLinearOperator)
-        $op(x, convert(AbstractMatrix, L))
-    end
-    @eval Base.$op(L::DiffEqArrayOperator, x::Number) = $op(convert(AbstractMatrix, L), x)
-    @eval Base.$op(x::Number, L::DiffEqArrayOperator) = $op(x, convert(AbstractMatrix, L))
-end
-
-### added in https://github.com/SciML/SciMLBase.jl/pull/377
-function LinearAlgebra.mul!(Y::AbstractVecOrMat, L::AbstractDiffEqLinearOperator,
-        B::AbstractVecOrMat)
-    mul!(Y, convert(AbstractMatrix, L), B)
-end
-###
-
-function LinearAlgebra.mul!(Y::AbstractArray, L::AbstractDiffEqLinearOperator,
-        B::AbstractArray)
-    mul!(Y, convert(AbstractMatrix, L), B)
-end
-
-### added in https://github.com/SciML/SciMLBase.jl/pull/377
-function LinearAlgebra.mul!(Y::AbstractVecOrMat, L::AbstractDiffEqLinearOperator,
-        B::AbstractVecOrMat, α::Number, β::Number)
-    mul!(Y, convert(AbstractMatrix, L), B, α, β)
-end
-###
-
-function LinearAlgebra.mul!(Y::AbstractArray, L::AbstractDiffEqLinearOperator,
-        B::AbstractArray, α::Number, β::Number)
-    mul!(Y, convert(AbstractMatrix, L), B, α, β)
-end
-
-for pred in (:isreal, :issymmetric, :ishermitian, :isposdef)
-    @eval function LinearAlgebra.$pred(L::AbstractDiffEqLinearOperator)
-        $pred(convert(AbstractArray, L))
-    end
-end
-for op in (:sum, :prod)
-    @eval function LinearAlgebra.$op(L::AbstractDiffEqLinearOperator; kwargs...)
-        $op(convert(AbstractArray, L); kwargs...)
-    end
-end
-function LinearAlgebra.factorize(L::AbstractDiffEqLinearOperator)
-    FactorizedDiffEqArrayOperator(factorize(convert(AbstractMatrix, L)))
-end
-for fact in (:lu, :lu!, :qr, :qr!, :cholesky, :cholesky!, :ldlt, :ldlt!,
-    :bunchkaufman, :bunchkaufman!, :lq, :lq!, :svd, :svd!)
-    @eval function LinearAlgebra.$fact(L::AbstractDiffEqLinearOperator, args...)
-        FactorizedDiffEqArrayOperator($fact(convert(AbstractMatrix, L), args...))
-    end
-    @eval function LinearAlgebra.$fact(L::AbstractDiffEqLinearOperator; kwargs...)
-        FactorizedDiffEqArrayOperator($fact(convert(AbstractMatrix, L); kwargs...))
-    end
-end
-
-# Routines that use the full matrix representation
-Base.Matrix(L::AbstractDiffEqLinearOperator) = Matrix(convert(AbstractMatrix, L))
-LinearAlgebra.exp(L::AbstractDiffEqLinearOperator) = exp(Matrix(L))
diff --git a/src/operators/diffeq_operator.jl b/src/operators/diffeq_operator.jl
deleted file mode 100644
index d31ebec3e7..0000000000
--- a/src/operators/diffeq_operator.jl
+++ /dev/null
@@ -1,180 +0,0 @@
-"""
-    AffineDiffEqOperator{T} <: AbstractDiffEqOperator{T}
-
-`Ex: (A₁(t) + ... + Aₙ(t))*u + B₁(t) + ... + Bₘ(t)`
-
-AffineDiffEqOperator{T}(As,Bs,du_cache=nothing)
-
-Takes in two tuples for split Affine DiffEqs
-
- 1. update_coefficients! works by updating the coefficients of the component
-    operators.
- 2. Function calls L(u, p, t) and L(du, u, p, t) are fallbacks interpreted in this form.
-    This will allow them to work directly in the nonlinear ODE solvers without
-    modification.
- 3. f(du, u, p, t) is only allowed if a du_cache is given
- 4. B(t) can be Union{Number,AbstractArray}, in which case they are constants.
-    Otherwise they are interpreted they are functions v=B(t) and B(v,t)
-
-Solvers will see this operator from integrator.f and can interpret it by
-checking the internals of As and Bs. For example, it can check isconstant(As[1])
-etc.
-"""
-struct AffineDiffEqOperator{T, T1, T2, U} <: AbstractDiffEqOperator{T}
-    As::T1
-    Bs::T2
-    du_cache::U
-    function AffineDiffEqOperator{T}(As, Bs, du_cache = nothing) where {T}
-        all([size(a) == size(As[1])
-             for a in As]) || error("Operator sizes do not agree")
-        new{T, typeof(As), typeof(Bs), typeof(du_cache)}(As, Bs, du_cache)
-    end
-end
-
-Base.size(L::AffineDiffEqOperator) = size(L.As[1])
-
-function (L::AffineDiffEqOperator)(u, p, t::Number)
-    update_coefficients!(L, u, p, t)
-    du = sum(A * u for A in L.As)
-    for B in L.Bs
-        if B isa Union{Number, AbstractArray}
-            du .+= B
-        else
-            du .+= B(t)
-        end
-    end
-    du
-end
-
-function (L::AffineDiffEqOperator)(du, u, p, t::Number)
-    update_coefficients!(L, u, p, t)
-    L.du_cache === nothing &&
-        error("Can only use inplace AffineDiffEqOperator if du_cache is given.")
-    du_cache = L.du_cache
-    fill!(du, zero(first(du)))
-    # TODO: Make type-stable via recursion
-    for A in L.As
-        mul!(du_cache, A, u)
-        du .+= du_cache
-    end
-    for B in L.Bs
-        if B isa Union{Number, AbstractArray}
-            du .+= B
-        else
-            B(du_cache, t)
-            du .+= du_cache
-        end
-    end
-end
-
-function update_coefficients!(L::AffineDiffEqOperator, u, p, t)
-    # TODO: Make type-stable via recursion
-    for A in L.As
-        update_coefficients!(A, u, p, t)
-    end
-    for B in L.Bs
-        update_coefficients!(B, u, p, t)
-    end
-end
-
-@deprecate is_constant(L::AbstractDiffEqOperator) isconstant(L)
-
-# Scaled operator (α * A)
-struct DiffEqScaledOperator{T, F, OpType <: AbstractDiffEqLinearOperator{T}} <:
-       AbstractDiffEqCompositeOperator{T}
-    coeff::DiffEqScalar{T, F}
-    op::OpType
-end
-
-# Recursive routines that use `getops`
-function update_coefficients!(L::AbstractDiffEqCompositeOperator, u, p, t)
-    for op in getops(L)
-        update_coefficients!(op, u, p, t)
-    end
-    L
-end
-
-function Base.:*(α::DiffEqScalar{T, F}, L::AbstractDiffEqLinearOperator{T}) where {T, F}
-    DiffEqScaledOperator(α, L)
-end
-function Base.:*(α::Number, L::AbstractDiffEqLinearOperator{T}) where {T}
-    DiffEqScaledOperator(DiffEqScalar(convert(T, α)), L)
-end
-Base.:-(L::AbstractDiffEqLinearOperator{T}) where {T} = DiffEqScalar(-one(T)) * L
-getops(L::DiffEqScaledOperator) = (L.coeff, L.op)
-Base.Matrix(L::DiffEqScaledOperator) = L.coeff * Matrix(L.op)
-function Base.convert(::Type{AbstractMatrix}, L::DiffEqScaledOperator)
-    L.coeff * convert(AbstractMatrix, L.op)
-end
-
-Base.size(L::DiffEqScaledOperator, i::Integer) = size(L.op, i)
-Base.size(L::DiffEqScaledOperator, args...) = size(L.op, args...)
-LinearAlgebra.opnorm(L::DiffEqScaledOperator, p::Real = 2) = abs(L.coeff) * opnorm(L.op, p)
-Base.getindex(L::DiffEqScaledOperator, i::Int) = L.coeff * L.op[i]
-Base.getindex(L::DiffEqScaledOperator, I::Vararg{Int, N}) where {N} = L.coeff * L.op[I...]
-
-Base.:*(L::DiffEqScaledOperator, x::AbstractVecOrMat) = L.coeff * (L.op * x)
-Base.:*(L::DiffEqScaledOperator, x::AbstractArray) = L.coeff * (L.op * x)
-
-Base.:*(x::AbstractVecOrMat, L::DiffEqScaledOperator) = (x * L.op) * L.coeff
-Base.:*(x::AbstractArray, L::DiffEqScaledOperator) = (x * L.op) * L.coeff
-
-function LinearAlgebra.mul!(r::AbstractVecOrMat, L::DiffEqScaledOperator,
-        x::AbstractVecOrMat)
-    mul!(r, L.op, x)
-    r .= r * L.coeff
-end
-function LinearAlgebra.mul!(r::AbstractArray, L::DiffEqScaledOperator, x::AbstractArray)
-    mul!(r, L.op, x)
-    r .= r * L.coeff
-end
-
-function LinearAlgebra.mul!(r::AbstractVecOrMat, x::AbstractVecOrMat,
-        L::DiffEqScaledOperator)
-    mul!(r, x, L.op)
-    r .= r * L.coeff
-end
-function LinearAlgebra.mul!(r::AbstractArray, x::AbstractArray, L::DiffEqScaledOperator)
-    mul!(r, x, L.op)
-    r .= r * L.coeff
-end
-
-Base.:/(L::DiffEqScaledOperator, x::AbstractVecOrMat) = L.coeff * (L.op / x)
-Base.:/(L::DiffEqScaledOperator, x::AbstractArray) = L.coeff * (L.op / x)
-
-Base.:/(x::AbstractVecOrMat, L::DiffEqScaledOperator) = 1 / L.coeff * (x / L.op)
-Base.:/(x::AbstractArray, L::DiffEqScaledOperator) = 1 / L.coeff * (x / L.op)
-
-Base.:\(L::DiffEqScaledOperator, x::AbstractVecOrMat) = 1 / L.coeff * (L.op \ x)
-Base.:\(L::DiffEqScaledOperator, x::AbstractArray) = 1 / L.coeff * (L.op \ x)
-
-Base.:\(x::AbstractVecOrMat, L::DiffEqScaledOperator) = L.coeff * (x \ L)
-Base.:\(x::AbstractArray, L::DiffEqScaledOperator) = L.coeff * (x \ L)
-
-for N in (2, 3)
-    @eval begin
-        function LinearAlgebra.mul!(Y::AbstractArray{T, $N},
-                L::DiffEqScaledOperator{T},
-                B::AbstractArray{T, $N}) where {T}
-            LinearAlgebra.lmul!(Y, L.coeff, mul!(Y, L.op, B))
-        end
-    end
-end
-
-function LinearAlgebra.ldiv!(Y::AbstractVecOrMat, L::DiffEqScaledOperator,
-        B::AbstractVecOrMat)
-    lmul!(1 / L.coeff, ldiv!(Y, L.op, B))
-end
-function LinearAlgebra.ldiv!(Y::AbstractArray, L::DiffEqScaledOperator, B::AbstractArray)
-    lmul!(1 / L.coeff, ldiv!(Y, L.op, B))
-end
-
-LinearAlgebra.factorize(L::DiffEqScaledOperator) = L.coeff * factorize(L.op)
-for fact in (:lu, :lu!, :qr, :qr!, :cholesky, :cholesky!, :ldlt, :ldlt!,
-    :bunchkaufman, :bunchkaufman!, :lq, :lq!, :svd, :svd!)
-    @eval function LinearAlgebra.$fact(L::DiffEqScaledOperator, args...)
-        L.coeff * fact(L.op, args...)
-    end
-end
-
-isconstant(L::AbstractDiffEqCompositeOperator) = all(isconstant, getops(L))
diff --git a/src/operators/operators.jl b/src/operators/operators.jl
deleted file mode 100644
index 6cd7774077..0000000000
--- a/src/operators/operators.jl
+++ /dev/null
@@ -1,15 +0,0 @@
-#
-# Extra standard assumptions
-isconstant(::AbstractDiffEqLinearOperator) = true
-islinear(o::AbstractDiffEqLinearOperator) = isconstant(o)
-
-# fallbacks
-islinear(L) = false
-isconstant(L) = false
-
-# Other ones from LinearMaps.jl
-# Generic fallbacks
-LinearAlgebra.exp(L::AbstractDiffEqLinearOperator, t) = exp(t * L)
-has_exp(L::AbstractDiffEqLinearOperator) = true
-expmv(L::AbstractDiffEqLinearOperator, u, p, t) = exp(L, t) * u
-expmv!(v, L::AbstractDiffEqLinearOperator, u, p, t) = mul!(v, exp(L, t), u)
diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl
index 2908ed75ad..463729570f 100644
--- a/src/scimlfunctions.jl
+++ b/src/scimlfunctions.jl
@@ -4220,7 +4220,7 @@ function BVPFunction{iip, specialize, twopoint}(f, bc;
         error("FunctionWrapperSpecialize must be used on the problem constructor for access to u0, p, and t types!")
     end
 
-    if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
+    if jac === nothing && isa(jac_prototype, AbstractSciMLOperator)
         if iip_f
             jac = update_coefficients! #(J,u,p,t)
         else
@@ -4228,7 +4228,7 @@ function BVPFunction{iip, specialize, twopoint}(f, bc;
         end
     end
 
-    if bcjac === nothing && isa(bcjac_prototype, AbstractDiffEqLinearOperator)
+    if bcjac === nothing && isa(bcjac_prototype, AbstractSciMLOperator)
         if iip_bc
             bcjac = update_coefficients! #(J,u,p,t)
         else
@@ -4393,7 +4393,7 @@ function DynamicalBVPFunction{iip, specialize, twopoint}(f, bc;
         error("FunctionWrapperSpecialize must be used on the problem constructor for access to u0, p, and t types!")
     end
 
-    if jac === nothing && isa(jac_prototype, AbstractDiffEqLinearOperator)
+    if jac === nothing && isa(jac_prototype, AbstractSciMLOperator)
         if iip_f
             jac = update_coefficients! #(J,u,p,t)
         else
@@ -4401,7 +4401,7 @@ function DynamicalBVPFunction{iip, specialize, twopoint}(f, bc;
         end
     end
 
-    if bcjac === nothing && isa(bcjac_prototype, AbstractDiffEqLinearOperator)
+    if bcjac === nothing && isa(bcjac_prototype, AbstractSciMLOperator)
         if iip_bc
             bcjac = update_coefficients! #(J,u,p,t)
         else
diff --git a/test/aqua.jl b/test/aqua.jl
index 2c7f9215c9..32708ce632 100644
--- a/test/aqua.jl
+++ b/test/aqua.jl
@@ -3,7 +3,7 @@ using SciMLBase
 using Aqua
 
 # https://github.com/JuliaArrays/FillArrays.jl/pull/163
-@test_broken isempty(detect_ambiguities(SciMLBase))
+@test isempty(detect_ambiguities(SciMLBase))
 
 @testset "Aqua tests (performance)" begin
     # This tests that we don't accidentally run into
diff --git a/test/diffeqoperator.jl b/test/diffeqoperator.jl
deleted file mode 100644
index bfd8d89bb2..0000000000
--- a/test/diffeqoperator.jl
+++ /dev/null
@@ -1,39 +0,0 @@
-using SciMLBase
-using LinearAlgebra
-
-@testset "DiffEqOperator" begin
-    n = 10
-    u = rand(n)
-    p = nothing
-    t = 0
-
-    A = rand(n, n)
-    At = A'
-
-    AA = SciMLBase.DiffEqArrayOperator(A)
-    AAt = AA'
-
-    @test AA isa SciMLBase.DiffEqArrayOperator
-    @test AAt isa SciMLBase.DiffEqArrayOperator
-
-    FF = factorize(AA)
-    FFt = FF'
-
-    @test FF isa SciMLBase.FactorizedDiffEqArrayOperator
-    @test FFt isa SciMLBase.FactorizedDiffEqArrayOperator
-
-    @test eachindex(A) === eachindex(AA)
-    @test eachindex(A') === eachindex(AAt) === eachindex(DiffEqArrayOperator(At))
-
-    @test A ≈ convert(AbstractMatrix, AA) ≈ convert(AbstractMatrix, FF)
-    @test At ≈ convert(AbstractMatrix, AAt) ≈ convert(AbstractMatrix, FFt)
-
-    @test A ≈ Matrix(AA) ≈ Matrix(FF)
-    @test At ≈ Matrix(AAt) ≈ Matrix(FFt)
-
-    @test A * u ≈ AA(u, p, t) ≈ FF(u, p, t)
-    @test At * u ≈ AAt(u, p, t) ≈ FFt(u, p, t)
-
-    @test A \ u ≈ AA \ u ≈ FF \ u
-    @test At \ u ≈ AAt \ u ≈ FFt \ u
-end
diff --git a/test/downstream/modelingtoolkit_remake.jl b/test/downstream/modelingtoolkit_remake.jl
index 5aefcaa3b6..ef1013b6b4 100644
--- a/test/downstream/modelingtoolkit_remake.jl
+++ b/test/downstream/modelingtoolkit_remake.jl
@@ -329,7 +329,7 @@ end
 
     newp = remake_buffer(sccprob.f.sys, sccprob.p, [σ], [3.0])
     sccprob4 = remake(sccprob; parameters_alias = false, p = newp,
-        probs = [remake(sccprob.probs[1]; p = [σ => 3.0]), sccprob.probs[2]])
+        probs = [remake(sccprob.probs[1]; p = deepcopy(newp)), sccprob.probs[2]])
     @test !sccprob4.parameters_alias
     @test sccprob4.p !== sccprob4.probs[1].p
     @test sccprob4.p !== sccprob4.probs[2].p
diff --git a/test/runtests.jl b/test/runtests.jl
index cba5acc9d5..766dbe6750 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -45,9 +45,6 @@ end
         @time @safetestset "Ensemble functionality" begin
             include("ensemble_tests.jl")
         end
-        @time @safetestset "DiffEqOperator tests" begin
-            include("diffeqoperator.jl")
-        end
         @time @safetestset "Solution interface" begin
             include("solution_interface.jl")
         end