diff --git a/Project.toml b/Project.toml
index a376f045..acf02c51 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ authors = ["Yi-Te Huang"]
 [deps]
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
@@ -30,6 +31,7 @@ HierarchicalEOM_CUDAExt = "CUDA"
 CUDA = "5"
 DiffEqCallbacks = "2 - 4"
 FastExpm = "1.1"
+FillArrays = "1"
 IncompleteLU = "0.2"
 LinearAlgebra = "1"
 LinearSolve = "2.4.2 - 3"
@@ -37,7 +39,7 @@ OrdinaryDiffEqCore = "1"
 OrdinaryDiffEqLowOrderRK = "1"
 Pkg = "1"
 ProgressMeter = "1.11.0"
-QuantumToolbox = "0.38.1"
+QuantumToolbox = "0.39"
 Reexport = "1"
 SciMLBase = "2.105"
 SciMLOperators = "1.4"
diff --git a/docs/make.jl b/docs/make.jl
index 0daab07f..03390e4c 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -6,7 +6,8 @@ using Documenter, HierarchicalEOM
 
 DocMeta.setdocmeta!(HierarchicalEOM, :DocTestSetup, :(using HierarchicalEOM); recursive = true)
 
-const DRAFT = false # set `true` to disable cell evaluation
+# some options for makedocs
+const DRAFT = get(ENV, "DRAFT", false) == "true"  # DRAFT   = true  disables cell evaluation
 
 const MathEngine = MathJax3(
     Dict(
diff --git a/docs/src/bath_boson/bosonic_bath_intro.md b/docs/src/bath_boson/bosonic_bath_intro.md
index a51373ff..f41a1c50 100644
--- a/docs/src/bath_boson/bosonic_bath_intro.md
+++ b/docs/src/bath_boson/bosonic_bath_intro.md
@@ -62,7 +62,7 @@ c_list = correlation_function(bath, tlist)
 ```
 Here, `c_list` is a list which contains the value of ``C(t)`` corresponds to the given time series `tlist`.
 
-## Exponent
+## Exponential terms of the Bosonic Bath
 `HierarchicalEOM.jl` also supports users to access the specific exponential term with brackets `[]`. This returns an [`Exponent`](@ref) object, which contains the corresponding value of ``\eta_i`` and ``\gamma_i``:
 ```julia
 e = bath[2] # the 2nd-term
diff --git a/docs/src/bath_boson_RWA/bosonic_bath_RWA_intro.md b/docs/src/bath_boson_RWA/bosonic_bath_RWA_intro.md
index 27dc993a..a9d3c4ec 100644
--- a/docs/src/bath_boson_RWA/bosonic_bath_RWA_intro.md
+++ b/docs/src/bath_boson_RWA/bosonic_bath_RWA_intro.md
@@ -47,7 +47,7 @@ cp_list, cm_list = correlation_function(bath, tlist)
 ```
 Here, `cp_list` and `cm_list` are the lists which contain the value of ``C^{\nu=+}(t)`` and ``C^{\nu=-}(t)`` correspond to the given time series `tlist`, respectively.
 
-## Exponent
+## Exponential terms of the Bosonic Bath (RWA)
 `HierarchicalEOM.jl` also supports users to access the specific exponential term with brackets `[]`. This returns an [`Exponent`](@ref) object, which contains the corresponding value of ``\eta_i^\nu`` and ``\gamma_i^\nu``:
 ```julia
 e = bath[2] # the 2nd-term
diff --git a/docs/src/bath_fermion/fermionic_bath_intro.md b/docs/src/bath_fermion/fermionic_bath_intro.md
index e85b87e7..6c8df65f 100644
--- a/docs/src/bath_fermion/fermionic_bath_intro.md
+++ b/docs/src/bath_fermion/fermionic_bath_intro.md
@@ -50,7 +50,7 @@ cp_list, cm_list = correlation_function(bath, tlist)
 ```
 Here, `cp_list` and `cm_list` are the lists which contain the value of ``C^{\nu=+}(t)`` and ``C^{\nu=-}(t)`` correspond to the given time series `tlist`, respectively.
 
-## Exponent
+## Exponential terms of the Fermionic Bath
 `HierarchicalEOM.jl` also supports users to access the specific exponential term with brackets `[]`. This returns an [`Exponent`](@ref) object, which contains the corresponding value of ``\eta_i^\nu`` and ``\gamma_i^\nu``:
 ```julia
 e = bath[2] # the 2nd-term
diff --git a/docs/src/libraryAPI.md b/docs/src/libraryAPI.md
index 5f652493..a74389bc 100644
--- a/docs/src/libraryAPI.md
+++ b/docs/src/libraryAPI.md
@@ -72,12 +72,14 @@ M_Fermion
 M_Fermion(Hsys::QuantumObject, tier::Int, Bath::Vector{FermionBath}, parity::AbstractParity=EVEN; threshold::Real=0.0, verbose::Bool=true)
 M_Boson_Fermion
 M_Boson_Fermion(Hsys::QuantumObject, tier_b::Int, tier_f::Int, Bath_b::Vector{BosonBath}, Bath_f::Vector{FermionBath}, parity::AbstractParity=EVEN; threshold::Real=0.0, verbose::Bool=true)
-size(M::HEOMSuperOp)
-size(M::HEOMSuperOp, dim::Int)
-size(M::AbstractHEOMLSMatrix)
-size(M::AbstractHEOMLSMatrix, dim::Int)
-eltype(M::HEOMSuperOp)
-eltype(M::AbstractHEOMLSMatrix)
+Base.size(M::HEOMSuperOp)
+Base.size(M::HEOMSuperOp, dim::Int)
+Base.size(M::AbstractHEOMLSMatrix)
+Base.size(M::AbstractHEOMLSMatrix, dim::Int)
+Base.eltype(M::HEOMSuperOp)
+Base.eltype(M::AbstractHEOMLSMatrix)
+iscached
+isconstant
 Propagator
 addBosonDissipator
 addFermionDissipator
@@ -88,8 +90,8 @@ addTerminator
 ```@docs
 ADOs
 ADOs(V::AbstractVector, N::Int)
-length(A::ADOs)
-eltype(A::ADOs)
+Base.length(A::ADOs)
+Base.eltype(A::ADOs)
 getRho
 getADO
 QuantumToolbox.expect
@@ -106,6 +108,7 @@ getIndexEnsemble
 ## [Time Evolution](@id lib-Time-Evolution)
 There are two function definitions of `HEOMsolve`, which depend on different methods to solve the time evolution:
 ```@docs
+HEOMsolveProblem
 HEOMsolve
 HEOMsolve_map
 TimeEvolutionHEOMSol
diff --git a/docs/src/time_evolution.md b/docs/src/time_evolution.md
index f704b10c..79f94fc0 100644
--- a/docs/src/time_evolution.md
+++ b/docs/src/time_evolution.md
@@ -66,7 +66,7 @@ The first method is implemented by solving the ordinary differential equation (O
 See the docstring of this method:  
 
 ```@docs
-HEOMsolve(M::AbstractHEOMLSMatrix, ρ0::T_state, tlist::AbstractVector; e_ops::Union{Nothing,AbstractVector} = nothing, alg::OrdinaryDiffEqAlgorithm = DP5(), H_t::Union{Nothing,Function} = nothing, params::NamedTuple = NamedTuple(), progress_bar::Union{Val,Bool} = Val(true), inplace::Union{Val,Bool} = Val(true), kwargs...) where {T_state<:Union{QuantumObject,ADOs}}
+HEOMsolve(M::AbstractHEOMLSMatrix, ρ0::T_state, tlist::AbstractVector; e_ops::Union{Nothing,AbstractVector} = nothing, alg::AbstractODEAlgorithm = DP5(), H_t::Union{Nothing,Function} = nothing, params::NamedTuple = NamedTuple(), progress_bar::Union{Val,Bool} = Val(true), inplace::Union{Val,Bool} = Val(true), kwargs...) where {T_state<:Union{QuantumObject,ADOs}}
 ```
 
 ```julia
diff --git a/src/HeomBase.jl b/src/HeomBase.jl
index e9b46290..261bdd19 100644
--- a/src/HeomBase.jl
+++ b/src/HeomBase.jl
@@ -1,5 +1,10 @@
 export AbstractHEOMLSMatrix
 
+@doc raw"""
+    abstract type AbstractHEOMLSMatrix{T}
+
+Abstract type for all HEOM liouvillian superoperators
+"""
 abstract type AbstractHEOMLSMatrix{T} end
 
 function Base.getproperty(M::AbstractHEOMLSMatrix, key::Symbol)
diff --git a/src/HierarchicalEOM.jl b/src/HierarchicalEOM.jl
index 595e8de7..7982aac9 100644
--- a/src/HierarchicalEOM.jl
+++ b/src/HierarchicalEOM.jl
@@ -4,6 +4,9 @@ module HierarchicalEOM
 using LinearAlgebra
 using SparseArrays
 
+import Base.Threads: @threads, nthreads, Channel
+import Pkg
+
 # Re-export QuantumToolbox
 import Reexport: @reexport
 @reexport using QuantumToolbox
@@ -17,11 +20,13 @@ import QuantumToolbox:
     _spost,
     _sprepost,
     _liouvillian,
+    _sum_lindblad_dissipators,
     _gen_dimensions,
     _get_dims_string,
     dimensions_to_dims,
     _save_func,
     _merge_saveat,
+    _merge_tstops,
     _merge_kwargs_with_callback,
     _get_expvals,
     _se_me_map_prob_func,
@@ -47,19 +52,24 @@ import SciMLBase:
     EnsembleThreads,
     FullSpecialize,
     CallbackSet,
-    NullParameters
+    NullParameters,
+    AbstractODEAlgorithm
 import SciMLOperators:
-    AbstractSciMLOperator, MatrixOperator, ScaledOperator, AddedOperator, update_coefficients!, concretize
-import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm
+    SciMLOperators,
+    AbstractSciMLOperator,
+    MatrixOperator,
+    ScaledOperator,
+    AddedOperator,
+    update_coefficients!,
+    concretize
 import OrdinaryDiffEqLowOrderRK: DP5
 import DiffEqCallbacks: FunctionCallingCallback, TerminateSteadyState
 import LinearSolve: LinearProblem, SciMLLinearSolveAlgorithm, KrylovJL_GMRES
 
 # other dependencies (in alphabetical order)
-import Base.Threads: @threads, nthreads, Channel
 import FastExpm: fastExpm
+import FillArrays: Eye
 import IncompleteLU: ilu
-import Pkg
 import ProgressMeter: Progress, next!
 
 # Basic functions
diff --git a/src/bath/BosonBath.jl b/src/bath/BosonBath.jl
index 7cb87e29..4399fd82 100644
--- a/src/bath/BosonBath.jl
+++ b/src/bath/BosonBath.jl
@@ -240,8 +240,7 @@ function bosonReal(op::QuantumObject, η_real::Vector{Ti}, γ_real::Vector{Tj})
     if N_exp_term != length(γ_real)
         error("The length of \'η_real\' and \'γ_real\' should be the same.")
     end
-    Id_cache = I(size(_op, 1))
-    return bosonReal(_spre(_op.data, Id_cache) - _spost(_op.data, Id_cache), _op.dimensions, η_real, γ_real, N_exp_term)
+    return bosonReal(_spre(_op.data) - _spost(_op.data), _op.dimensions, η_real, γ_real, N_exp_term)
 end
 
 @doc raw"""
@@ -285,9 +284,8 @@ function bosonImag(op::QuantumObject, η_imag::Vector{Ti}, γ_imag::Vector{Tj})
     if N_exp_term != length(γ_imag)
         error("The length of \'η_imag\' and \'γ_imag\' should be the same.")
     end
-    Id_cache = I(size(_op, 1))
-    spreQ = _spre(_op.data, Id_cache)
-    spostQ = _spost(_op.data, Id_cache)
+    spreQ = _spre(_op.data)
+    spostQ = _spost(_op.data)
     return bosonImag(spreQ - spostQ, spreQ + spostQ, _op.dimensions, η_imag, γ_imag, N_exp_term)
 end
 
@@ -340,9 +338,8 @@ function bosonRealImag(
     if (N_exp_term != length(η_imag)) || (N_exp_term != length(γ))
         error("The length of \'η_real\', \'η_imag\' and \'γ\' should be the same.")
     end
-    Id_cache = I(size(_op, 1))
-    spreQ = _spre(_op.data, Id_cache)
-    spostQ = _spost(_op.data, Id_cache)
+    spreQ = _spre(_op.data)
+    spostQ = _spost(_op.data)
     return bosonRealImag(spreQ - spostQ, spreQ + spostQ, _op.dimensions, η_real, η_imag, γ, N_exp_term)
 end
 
@@ -437,11 +434,10 @@ function bosonAbsorb(
     if (N_exp_term != length(γ_absorb)) || (N_exp_term != length(η_emit))
         error("The length of \'η_absorb\', \'γ_absorb\' and \'η_emit\' should all be the same.")
     end
-    Id_cache = I(size(_op, 1))
     return bosonAbsorb(
-        _spre(_op.data, Id_cache),
-        _spost(_op.data, Id_cache),
-        _spre(adjoint(_op).data, Id_cache) - _spost(adjoint(_op).data, Id_cache),
+        _spre(_op.data),
+        _spost(_op.data),
+        _spre(adjoint(_op).data) - _spost(adjoint(_op).data),
         _op.dimensions,
         η_absorb,
         γ_absorb,
@@ -501,11 +497,10 @@ function bosonEmit(
     if (N_exp_term != length(γ_emit)) || (N_exp_term != length(η_absorb))
         error("The length of \'η_emit\', \'γ_emit\' and \'η_absorb\' should all be the same.")
     end
-    Id_cache = I(size(_op, 1))
     return bosonEmit(
-        _spre(_op.data, Id_cache),
-        _spost(_op.data, Id_cache),
-        _spre(adjoint(_op).data, Id_cache) - _spost(adjoint(_op).data, Id_cache),
+        _spre(_op.data),
+        _spost(_op.data),
+        _spre(adjoint(_op).data) - _spost(adjoint(_op).data),
         _op.dimensions,
         η_emit,
         γ_emit,
diff --git a/src/bath/FermionBath.jl b/src/bath/FermionBath.jl
index 738e28d0..056ce892 100644
--- a/src/bath/FermionBath.jl
+++ b/src/bath/FermionBath.jl
@@ -136,12 +136,11 @@ function fermionAbsorb(
     if (N_exp_term != length(γ_absorb)) || (N_exp_term != length(η_emit))
         error("The length of \'η_absorb\', \'γ_absorb\' and \'η_emit\' should all be the same.")
     end
-    Id_cache = I(size(_op, 1))
     return fermionAbsorb(
-        _spre(_op.data, Id_cache),
-        _spost(_op.data, Id_cache),
-        _spre(adjoint(_op).data, Id_cache),
-        _spost(adjoint(_op).data, Id_cache),
+        _spre(_op.data),
+        _spost(_op.data),
+        _spre(adjoint(_op).data),
+        _spost(adjoint(_op).data),
         _op.dimensions,
         η_absorb,
         γ_absorb,
@@ -204,12 +203,11 @@ function fermionEmit(
         error("The length of \'η_emit\', \'γ_emit\' and \'η_absorb\' should all be the same.")
     end
 
-    Id_cache = I(size(_op, 1))
     return fermionEmit(
-        _spre(_op.data, Id_cache),
-        _spost(_op.data, Id_cache),
-        _spre(adjoint(_op).data, Id_cache),
-        _spost(adjoint(_op).data, Id_cache),
+        _spre(_op.data),
+        _spost(_op.data),
+        _spre(adjoint(_op).data),
+        _spost(adjoint(_op).data),
         _op.dimensions,
         η_emit,
         γ_emit,
diff --git a/src/density_of_states.jl b/src/density_of_states.jl
index 7ca724e2..ab3201f8 100644
--- a/src/density_of_states.jl
+++ b/src/density_of_states.jl
@@ -53,10 +53,8 @@ Calculate density of states for the fermionic system in frequency domain.
 
     # Handle d_op
     _tr = _Tr(M)
-    Id_sys = I(prod(d_op.dimensions))
-    Id_HEOM = I(M.N)
-    d_normal = HEOMSuperOp(spre(d_op, Id_sys), ODD, M; Id_cache = Id_HEOM)
-    d_dagger = HEOMSuperOp(spre(d_op', Id_sys), ODD, M; Id_cache = Id_HEOM)
+    d_normal = HEOMSuperOp(spre(d_op), ODD, M)
+    d_dagger = HEOMSuperOp(spre(d_op'), ODD, M)
     b_m = _HandleVectorType(M, (d_normal * ados).data)
     b_p = _HandleVectorType(M, (d_dagger * ados).data)
     _tr_d_normal = _HandleTraceVectorType(M, adjoint(d_normal.data) * _tr) # another adjoint will be applied in dot function later
@@ -80,7 +78,7 @@ Calculate density of states for the fermionic system in frequency domain.
         QuantumToolbox.settings.ProgressMeterKWARGS...,
     )
     i = convert(ElType, 1im)
-    I_total = I(size(M, 1))
+    I_total = Eye(size(M, 1))
     Iω1 = i * ωList[1] * I_total
     cache_m = init(LinearProblem(M.data.A - Iω1, b_m), alg, kwargs...)
     cache_p = init(LinearProblem(M.data.A + Iω1, b_p), alg, kwargs...)
diff --git a/src/deprecated.jl b/src/deprecated.jl
index b02f3645..b23efcd4 100644
--- a/src/deprecated.jl
+++ b/src/deprecated.jl
@@ -56,14 +56,9 @@ _HEOMSuperOp_deprecated_method_error() = error(
     "Specifying `mul_basis = \"L\", \"R\", or \"LR\"` to `HEOMSuperOp` is deprecated, try to construct system SuperOperator manually by using `spre`, `spost`, or `sprepost`.",
 )
 
-HEOMSuperOp(op, opParity::AbstractParity, dims, N::Int, mul_basis::AbstractString; Id_cache = I(N)) =
+HEOMSuperOp(op, opParity::AbstractParity, dims, N::Int, mul_basis::AbstractString) =
     _HEOMSuperOp_deprecated_method_error()
-HEOMSuperOp(
-    op,
-    opParity::AbstractParity,
-    refHEOMLS::AbstractHEOMLSMatrix,
-    mul_basis::AbstractString;
-    Id_cache = I(refHEOMLS.N),
-) = HEOMSuperOp(op, opParity, refHEOMLS.dimensions, refHEOMLS.N, mul_basis; Id_cache = Id_cache)
-HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs, mul_basis::AbstractString; Id_cache = I(refADOs.N)) =
-    HEOMSuperOp(op, opParity, refADOs.dimensions, refADOs.N, mul_basis; Id_cache = Id_cache)
+HEOMSuperOp(op, opParity::AbstractParity, refHEOMLS::AbstractHEOMLSMatrix, mul_basis::AbstractString) =
+    HEOMSuperOp(op, opParity, refHEOMLS.dimensions, refHEOMLS.N, mul_basis)
+HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs, mul_basis::AbstractString) =
+    HEOMSuperOp(op, opParity, refADOs.dimensions, refADOs.N, mul_basis)
diff --git a/src/evolution.jl b/src/evolution.jl
index 5ed25ccf..140a042b 100644
--- a/src/evolution.jl
+++ b/src/evolution.jl
@@ -1,4 +1,4 @@
-export HEOMsolve, heomsolve, HEOMsolve_map, heomsolve_map, TimeEvolutionHEOMSol
+export HEOMsolve, heomsolve, HEOMsolve_map, heomsolve_map, TimeEvolutionHEOMSol, HEOMsolveProblem
 
 @doc raw"""
     struct TimeEvolutionHEOMSol
@@ -23,7 +23,7 @@ struct TimeEvolutionHEOMSol{
     TS<:Vector{ADOs},
     TE<:Union{Nothing,AbstractMatrix},
     RETT<:Union{Nothing,Enum},
-    AlgT<:Union{Nothing,OrdinaryDiffEqAlgorithm},
+    AlgT<:Union{Nothing,AbstractODEAlgorithm},
     TolT<:Union{Nothing,Real},
 }
     Btier::Int
@@ -104,10 +104,8 @@ function HEOMsolve(
         expvals = Array{ComplexF64}(undef, 0, steps + 1)
         is_empty_e_ops = true
     else
-        Id_sys = I(prod(M.dimensions))
-        Id_HEOM = I(M.N)
         expvals = Array{ComplexF64}(undef, length(e_ops), steps + 1)
-        tr_e_ops = _generate_Eops(M, e_ops, Id_sys, Id_HEOM)
+        tr_e_ops = _generate_Eops(M, e_ops)
         is_empty_e_ops = isempty(e_ops)
     end
 
@@ -172,6 +170,44 @@ function _gen_ados_ode_vector(ados::ADOs, M::AbstractHEOMLSMatrix)
     return _HandleVectorType(M, ados.data)
 end
 
+@doc raw"""
+    HEOMsolveProblem(
+        M::AbstractHEOMLSMatrix,
+        ρ0::Union{QuantumObject,ADOs},
+        tlist::AbstractVector;
+        e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
+        H_t::Union{Nothing,QuantumObjectEvolution} = nothing,
+        params = NullParameters(),
+        progress_bar::Union{Val,Bool} = Val(true),
+        inplace::Union{Val,Bool} = Val(true),
+        kwargs...,
+    )
+
+Generate the ODEProblem for the time evolution of auxiliary density operators.
+
+# Parameters
+- `M::AbstractHEOMLSMatrix` : the matrix given from HEOM model
+- `ρ0::Union{QuantumObject,ADOs}` : system initial state (density matrix) or initial auxiliary density operators (`ADOs`)
+- `tlist::AbstractVector` : Denote the specific time points to save the solution at, during the solving process.
+- `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
+- `alg::AbstractODEAlgorithm` : The solving algorithm in package `DifferentialEquations.jl`. Default to `DP5()`.
+- `H_t::Union{Nothing,QuantumObjectEvolution}`: The time-dependent system Hamiltonian or Liouvillian. Default to `nothing`.
+- `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
+- `progress_bar::Union{Val,Bool}`: Whether to show the progress bar. Defaults to `Val(true)`. Using non-`Val` types might lead to type instabilities.
+- `inplace::Union{Val,Bool}`: Whether to use the inplace version of the ODEProblem. Defaults to `Val(true)`.
+- `kwargs` : The keyword arguments for the `ODEProblem`.
+
+# Notes
+- The [`ADOs`](@ref) will be saved depend on the keyword argument `saveat` in `kwargs`.
+- If `e_ops` is specified, the default value of `saveat=[tlist[end]]` (only save the final `ADOs`), otherwise, `saveat=tlist` (saving the `ADOs` corresponding to `tlist`). You can also specify `e_ops` and `saveat` separately.
+- The default tolerances in `kwargs` are given as `reltol=1e-6` and `abstol=1e-8`.
+- For more details about `alg` please refer to [`DifferentialEquations.jl` (ODE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
+- For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/)
+
+# Returns
+
+- `prob`: The `TimeEvolutionProblem` containing the `ODEProblem` for the time evolution of auxiliary density operators.
+"""
 function HEOMsolveProblem(
     M::AbstractHEOMLSMatrix,
     ρ0::T_state,
@@ -192,12 +228,12 @@ function HEOMsolveProblem(
     # handle initial state
     u0 = _gen_ados_ode_vector(ρ0, M)
 
-    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...)
-    kwargs3 = _generate_heom_kwargs(e_ops, makeVal(progress_bar), tlist, kwargs2, SaveFuncHEOMSolve, M)
-
     # define ODE problem (L should be an AbstractSciMLOperator)
     L = _make_L(M, H_t)
-    prob = ODEProblem{getVal(inplace),FullSpecialize}(L, u0, tspan, params; kwargs3...)
+    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...)
+    kwargs3 = _merge_tstops(kwargs2, isconstant(L), tlist)
+    kwargs4 = _generate_heom_kwargs(e_ops, makeVal(progress_bar), tlist, kwargs3, SaveFuncHEOMSolve, M)
+    prob = ODEProblem{getVal(inplace),FullSpecialize}(L, u0, tspan, params; kwargs4...)
 
     return TimeEvolutionProblem(prob, tlist, M.dimensions, (M = M,))
 end
@@ -213,7 +249,7 @@ Solve the time evolution for auxiliary density operators based on ordinary diffe
 - `ρ0::Union{QuantumObject,ADOs}` : system initial state (density matrix) or initial auxiliary density operators (`ADOs`)
 - `tlist::AbstractVector` : Denote the specific time points to save the solution at, during the solving process.
 - `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
-- `alg::OrdinaryDiffEqAlgorithm` : The solving algorithm in package `DifferentialEquations.jl`. Default to `DP5()`.
+- `alg::AbstractODEAlgorithm` : The solving algorithm in package `DifferentialEquations.jl`. Default to `DP5()`.
 - `H_t::Union{Nothing,QuantumObjectEvolution}`: The time-dependent system Hamiltonian or Liouvillian. Default to `nothing`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
 - `progress_bar::Union{Val,Bool}`: Whether to show the progress bar. Defaults to `Val(true)`. Using non-`Val` types might lead to type instabilities.
@@ -238,7 +274,7 @@ function HEOMsolve(
     ρ0::T_state,
     tlist::AbstractVector;
     e_ops::Union{Nothing,AbstractVector} = nothing,
-    alg::OrdinaryDiffEqAlgorithm = DP5(),
+    alg::AbstractODEAlgorithm = DP5(),
     H_t::Union{Nothing,QuantumObjectEvolution} = nothing,
     params = NullParameters(),
     progress_bar::Union{Val,Bool} = Val(true),
@@ -275,7 +311,7 @@ function HEOMsolve(
     end
 end
 
-function HEOMsolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = DP5(); kwargs...)
+function HEOMsolve(prob::TimeEvolutionProblem, alg::AbstractODEAlgorithm = DP5(); kwargs...)
     sol = solve(prob.prob, alg; kwargs...)
 
     return _gen_HEOMsolve_solution(sol, prob.times, prob.kwargs.M)
@@ -302,13 +338,10 @@ end
 
 const heomsolve = HEOMsolve # a synonym to align with qutip
 
-function _generate_Eops(M::AbstractHEOMLSMatrix, e_ops, Id_sys, Id_HEOM)
+function _generate_Eops(M::AbstractHEOMLSMatrix, e_ops)
     tr_e_ops = [
         # another adjoint will be applied in dot function in the HEOMsolveCallback
-        _HandleTraceVectorType(
-            M,
-            adjoint(HEOMSuperOp(spre(op, Id_sys), EVEN, M.dimensions, M.N; Id_cache = Id_HEOM).data) * _Tr(M),
-        ) for op in e_ops
+        _HandleTraceVectorType(M, adjoint(HEOMSuperOp(spre(op), EVEN, M.dimensions, M.N).data) * _Tr(M)) for op in e_ops
     ]
     return tr_e_ops
 end
@@ -338,7 +371,7 @@ function _generate_heom_kwargs(
     method::Type{SaveFuncHEOMSolve},
     M::AbstractHEOMLSMatrix,
 )
-    tr_e_ops = e_ops isa Nothing ? nothing : _generate_Eops(M, e_ops, I(prod(M.dimensions)), I(M.N))
+    tr_e_ops = e_ops isa Nothing ? nothing : _generate_Eops(M, e_ops)
 
     progr =
         getVal(progress_bar) ?
@@ -367,19 +400,18 @@ _generate_heom_kwargs(
 
 _make_L(M::AbstractHEOMLSMatrix, H_t::Nothing) = M.data
 function _make_L(M::AbstractHEOMLSMatrix, H_t::QuantumObjectEvolution)
-    Id_HEOM = I(M.N)
     MType = _get_SciML_matrix_wrapper(M)
     L_t = HandleMatrixType(liouvillian(H_t), M.dimensions, "H_t"; type = SuperOperator())
 
-    return M.data + _L_t_to_HEOMSuperOp(MType, L_t.data, Id_HEOM)
+    return M.data + _L_t_to_HEOMSuperOp(MType, L_t.data, M.N)
 end
 
-_L_t_to_HEOMSuperOp(MType::Type{<:AbstractSparseMatrix}, M::MatrixOperator, Id::Diagonal) =
-    MatrixOperator(MType(kron(Id, M.A)))
-_L_t_to_HEOMSuperOp(MType::Type{<:AbstractSparseMatrix}, M::ScaledOperator, Id::Diagonal) =
-    ScaledOperator(M.λ, _L_t_to_HEOMSuperOp(MType, M.L, Id))
-_L_t_to_HEOMSuperOp(MType::Type{<:AbstractSparseMatrix}, M::AddedOperator, Id::Diagonal) =
-    AddedOperator(map(op -> _L_t_to_HEOMSuperOp(MType, op, Id), M.ops))
+_L_t_to_HEOMSuperOp(MType::Type{<:AbstractSparseMatrix}, M::MatrixOperator, N::Int) =
+    MatrixOperator(MType(kron(Eye(N), M.A)))
+_L_t_to_HEOMSuperOp(MType::Type{<:AbstractSparseMatrix}, M::ScaledOperator, N::Int) =
+    ScaledOperator(M.λ, _L_t_to_HEOMSuperOp(MType, M.L, N))
+_L_t_to_HEOMSuperOp(MType::Type{<:AbstractSparseMatrix}, M::AddedOperator, N::Int) =
+    AddedOperator(map(op -> _L_t_to_HEOMSuperOp(MType, op, N), M.ops))
 
 @doc raw"""
     HEOMsolve_map(M, ρ0, tlist; e_ops, alg, ensemblealg, H_t, params, progress_bar, kwargs...)
@@ -394,7 +426,7 @@ This function computes the time evolution for all combinations (Cartesian produc
 - `ρ0::AbstractVector{<:Union{QuantumObject,ADOs}}` : system initial state(s) [the elements can be density matrix or `ADOs`].
 - `tlist::AbstractVector` : Denote the specific time points to save the solution at, during the solving process.
 - `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
-- `alg::OrdinaryDiffEqAlgorithm` : The solving algorithm in package `DifferentialEquations.jl`. Default to `DP5()`.
+- `alg::AbstractODEAlgorithm` : The solving algorithm in package `DifferentialEquations.jl`. Default to `DP5()`.
 - `ensemblealg::EnsembleAlgorithm`: Ensemble algorithm to use for parallel computation. Default is `EnsembleThreads()`.
 - `H_t::Union{Nothing,QuantumObjectEvolution}`: The time-dependent system Hamiltonian or Liouvillian. Default to `nothing`.
 - `params::Union{NullParameters,Tuple}`: A `Tuple` of parameter sets. Each element should be an `AbstractVector` representing the sweep range for that parameter. The function will solve for all combinations of initial states and parameter sets.
@@ -417,7 +449,7 @@ function HEOMsolve_map(
     ρ0::AbstractVector{<:T_state},
     tlist::AbstractVector;
     e_ops::Union{Nothing,AbstractVector} = nothing,
-    alg::OrdinaryDiffEqAlgorithm = DP5(),
+    alg::AbstractODEAlgorithm = DP5(),
     ensemblealg::EnsembleAlgorithm = EnsembleThreads(),
     H_t::Union{Nothing,QuantumObjectEvolution} = nothing,
     params::Union{NullParameters,Tuple} = NullParameters(),
@@ -463,7 +495,7 @@ HEOMsolve_map(
 function HEOMsolve_map(
     prob::TimeEvolutionProblem{<:ODEProblem},
     iter::AbstractArray,
-    alg::OrdinaryDiffEqAlgorithm = DP5(),
+    alg::AbstractODEAlgorithm = DP5(),
     ensemblealg::EnsembleAlgorithm = EnsembleThreads();
     prob_func::Union{Function,Nothing} = nothing,
     output_func::Union{Tuple,Nothing} = nothing,
diff --git a/src/heom_matrices/heom_matrix_base.jl b/src/heom_matrices/heom_matrix_base.jl
index 0aceffdb..def58c7f 100644
--- a/src/heom_matrices/heom_matrix_base.jl
+++ b/src/heom_matrices/heom_matrix_base.jl
@@ -32,7 +32,7 @@ function Base.getproperty(M::HEOMSuperOp, key::Symbol)
 end
 
 @doc raw"""
-    HEOMSuperOp(op, opParity, refHEOMLS; Id_cache=I(refHEOMLS.N))
+    HEOMSuperOp(op, opParity, refHEOMLS)
 Construct the HEOM superoperator matrix corresponding to the given system SuperOperator which acts on all `ADOs`.  
 
 During the multiplication on all the `ADOs`, the parity of the output `ADOs` might change depend on the parity of this HEOM superoperator.
@@ -42,11 +42,11 @@ During the multiplication on all the `ADOs`, the parity of the output `ADOs` mig
 - `opParity::AbstractParity` : the parity label of the given operator (`op`), should be `EVEN` or `ODD`.
 - `refHEOMLS::AbstractHEOMLSMatrix` : copy the system `dimensions` and number of `ADOs` (`N`) from this reference HEOMLS matrix
 """
-HEOMSuperOp(op, opParity::AbstractParity, refHEOMLS::AbstractHEOMLSMatrix; Id_cache = I(refHEOMLS.N)) =
-    HEOMSuperOp(op, opParity, refHEOMLS.dimensions, refHEOMLS.N; Id_cache = Id_cache)
+HEOMSuperOp(op, opParity::AbstractParity, refHEOMLS::AbstractHEOMLSMatrix) =
+    HEOMSuperOp(op, opParity, refHEOMLS.dimensions, refHEOMLS.N)
 
 @doc raw"""
-    HEOMSuperOp(op, opParity, refADOs; Id_cache=I(refADOs.N))
+    HEOMSuperOp(op, opParity, refADOs)
 Construct the HEOM superoperator matrix corresponding to the given system SuperOperator which acts on all `ADOs`.  
 
 During the multiplication on all the `ADOs`, the parity of the output `ADOs` might change depend on the parity of this HEOM superoperator.
@@ -56,11 +56,10 @@ During the multiplication on all the `ADOs`, the parity of the output `ADOs` mig
 - `opParity::AbstractParity` : the parity label of the given operator (`op`), should be `EVEN` or `ODD`.
 - `refADOs::ADOs` : copy the system `dimensions` and number of `ADOs` (`N`) from this reference `ADOs`   
 """
-HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs; Id_cache = I(refADOs.N)) =
-    HEOMSuperOp(op, opParity, refADOs.dimensions, refADOs.N; Id_cache = Id_cache)
+HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs) = HEOMSuperOp(op, opParity, refADOs.dimensions, refADOs.N)
 
 @doc raw"""
-    HEOMSuperOp(op, opParity, dims, N; Id_cache=I(N))
+    HEOMSuperOp(op, opParity, dims, N)
 Construct the HEOM superoperator matrix corresponding to the given system SuperOperator which acts on all `ADOs`.  
 
 During the multiplication on all the `ADOs`, the parity of the output `ADOs` might change depend on the parity of this HEOM superoperator.
@@ -71,11 +70,11 @@ During the multiplication on all the `ADOs`, the parity of the output `ADOs` mig
 - `dims` : the dimension list of the coupling operator (should be equal to the system `dimensions`).
 - `N::Int` : the number of `ADOs`.
 """
-function HEOMSuperOp(op, opParity::AbstractParity, dims, N::Int; Id_cache = I(N))
+function HEOMSuperOp(op, opParity::AbstractParity, dims, N::Int)
     dimensions = _gen_dimensions(dims)
     sup_op = HandleMatrixType(op, dimensions, "op (operator)"; type = SuperOperator())
 
-    return HEOMSuperOp(kron(Id_cache, sup_op.data), dimensions, N, opParity)
+    return HEOMSuperOp(kron(Eye(N), sup_op.data), dimensions, N, opParity)
 end
 
 @doc raw"""
@@ -203,6 +202,20 @@ function Base.:(-)(M::AbstractHEOMLSMatrix, Sup::HEOMSuperOp)
     return _reset_HEOMLS_data(M, M.data - Sup)
 end
 
+@doc raw"""
+    SciMLOperators.iscached(M::AbstractHEOMLSMatrix)
+
+Test whether the [`AbstractHEOMLSMatrix`](@ref) `M` has preallocated caches for inplace evaluations.
+"""
+SciMLOperators.iscached(M::AbstractHEOMLSMatrix) = iscached(M.data)
+
+@doc raw"""
+    SciMLOperators.isconstant(M::AbstractHEOMLSMatrix)
+
+Test whether the [`AbstractHEOMLSMatrix`](@ref) `M` is constant in time.
+"""
+SciMLOperators.isconstant(M::AbstractHEOMLSMatrix) = isconstant(M.data)
+
 @doc raw"""
     Propagator(M, Δt; threshold, nonzero_tol)
 Use `FastExpm.jl` to calculate the propagator matrix from a given HEOM Liouvillian superoperator matrix ``M`` with a specific time step ``\Delta t``.
@@ -271,10 +284,7 @@ Note that if ``V`` is acting on fermionic systems, it should be even-parity to b
 """
 function addBosonDissipator(M::AbstractHEOMLSMatrix, jumpOP::Vector{T} = QuantumObject[]) where {T<:QuantumObject}
     if length(jumpOP) > 0
-        Id_cache = I(prod(M.dimensions))
-        L = mapreduce(J -> lindblad_dissipator(J, Id_cache), +, jumpOP)
-
-        return M + HEOMSuperOp(L, M.parity, M)
+        return M + HEOMSuperOp(_sum_lindblad_dissipators(jumpOP), M.parity, M)
     else
         return M
     end
@@ -307,8 +317,7 @@ Note that the parity of the dissipator will be determined by the parity of the g
 """
 function addFermionDissipator(M::AbstractHEOMLSMatrix, jumpOP::Vector{T} = QuantumObject[]) where {T<:QuantumObject}
     if length(jumpOP) > 0
-        Id_cache = I(prod(M.dimensions))
-        L_data = mapreduce(J -> _fermion_lindblad_dissipator(J, M.parity, Id_cache), +, jumpOP)
+        L_data = mapreduce(J -> _fermion_lindblad_dissipator(J, M.parity), +, jumpOP)
         L = QuantumObject(L_data, type = SuperOperator(), dims = M.dimensions)
 
         return M + HEOMSuperOp(L, M.parity, M)
@@ -318,10 +327,10 @@ function addFermionDissipator(M::AbstractHEOMLSMatrix, jumpOP::Vector{T} = Quant
 end
 addFermionDissipator(M::AbstractHEOMLSMatrix, jumpOP::QuantumObject) = addFermionDissipator(M, [jumpOP])
 
-function _fermion_lindblad_dissipator(J::QuantumObject{Operator}, parity::AbstractParity, Id_cache = I(size(J, 1)))
+function _fermion_lindblad_dissipator(J::QuantumObject{Operator}, parity::AbstractParity)
     _J = J.data
     Jd_J = _J' * _J
-    return (-1)^(value(parity)) * _sprepost(_J, _J') - (_spre(Jd_J, Id_cache) + _spost(Jd_J, Id_cache)) / 2
+    return (-1)^(value(parity)) * _sprepost(_J, _J') - (_spre(Jd_J) + _spost(Jd_J)) / 2
 end
 
 @doc raw"""
@@ -435,7 +444,7 @@ function bath_sum_γ(nvec, baths::Vector{T}) where {T<:Union{AbstractBosonBath,A
 end
 
 # commutator of system Hamiltonian
-minus_i_L_op(Hsys::QuantumObject, Id = I(size(Hsys, 1))) = liouvillian(Hsys, Id).data
+minus_i_L_op(Hsys::QuantumObject) = liouvillian(Hsys).data
 
 # connect to bosonic (n-1)th-level for "Real & Imag combined operator"
 minus_i_D_op(bath::bosonRealImag, k, n_k) = n_k * (-1.0im * bath.η_real[k] * bath.Comm + bath.η_imag[k] * bath.anComm)
diff --git a/src/power_spectrum.jl b/src/power_spectrum.jl
index 5fd4ab88..3f741490 100644
--- a/src/power_spectrum.jl
+++ b/src/power_spectrum.jl
@@ -89,14 +89,12 @@ remember to set the parameters:
     end
     _check_sys_dim_and_ADOs_num(M, ados)
 
-    Id_HEOM = I(M.N)
-
     # Handle P_op
     if P_op isa HEOMSuperOp
         _check_sys_dim_and_ADOs_num(M, P_op)
         _P = P_op
     else
-        _P = HEOMSuperOp(spre(P_op), EVEN, M; Id_cache = Id_HEOM)
+        _P = HEOMSuperOp(spre(P_op), EVEN, M)
     end
     _tr_P = _HandleTraceVectorType(M, adjoint(_P.data) * _Tr(M)) # another adjoint will be applied in dot function later
 
@@ -111,7 +109,7 @@ remember to set the parameters:
         else
             new_parity = !ados.parity
         end
-        _Q_ados = HEOMSuperOp(spre(Q_op), new_parity, M; Id_cache = Id_HEOM) * ados
+        _Q_ados = HEOMSuperOp(spre(Q_op), new_parity, M) * ados
     end
     b = _HandleVectorType(M, _Q_ados.data)
 
@@ -133,7 +131,7 @@ remember to set the parameters:
         QuantumToolbox.settings.ProgressMeterKWARGS...,
     )
     i = reverse ? convert(ElType, 1im) : i = convert(ElType, -1im)
-    I_total = I(size(M, 1))
+    I_total = Eye(size(M, 1))
     cache = init(LinearProblem(M.data.A + i * ωList[1] * I_total, b), alg, kwargs...)
     for (idx, ω) in enumerate(ωList)
         if idx > 1
diff --git a/src/steadystate.jl b/src/steadystate.jl
index ff6c54c0..3c16c716 100644
--- a/src/steadystate.jl
+++ b/src/steadystate.jl
@@ -70,7 +70,7 @@ Solve the steady state of the auxiliary density operators based on time evolutio
 - `M::AbstractHEOMLSMatrix` : the matrix given from HEOM model, where the parity should be `EVEN`.
 - `ρ0::Union{QuantumObject,ADOs}` : system initial state (density matrix) or initial auxiliary density operators (`ADOs`)
 - `tspan::Number` : the time limit to find stationary state. Default to `Inf`
-- `alg::OrdinaryDiffEqAlgorithm` : The ODE algorithm in package `DifferentialEquations.jl`. Default to `DP5()`.
+- `alg::AbstractODEAlgorithm` : The ODE algorithm in package `DifferentialEquations.jl`. Default to `DP5()`.
 - `verbose::Bool` : To display verbose output or not. Defaults to `true`.
 - `kwargs` : The keyword arguments in `ODEProblem`
 
@@ -85,7 +85,7 @@ function QuantumToolbox.steadystate(
     M::AbstractHEOMLSMatrix{<:MatrixOperator},
     ρ0::T_state,
     tspan::Number = Inf;
-    alg::OrdinaryDiffEqAlgorithm = DP5(),
+    alg::AbstractODEAlgorithm = DP5(),
     verbose::Bool = true,
     kwargs...,
 ) where {T_state<:Union{QuantumObject,ADOs}}
diff --git a/test/M_Boson.jl b/test/M_Boson.jl
index a7353a14..4a78267a 100644
--- a/test/M_Boson.jl
+++ b/test/M_Boson.jl
@@ -32,6 +32,8 @@
     @test nnz(L.data.A) == nnz(L(0)) == 4422
     L = addBosonDissipator(L, J)
     @test nnz(L.data.A) == nnz(L(0)) == 4760
+    @test isconstant(L)
+    @test iscached(L)
     ados = steadystate(L; verbose = false)
     @test ados.dims == L.dims
     @test length(ados) == L.N
diff --git a/test/M_Boson_Fermion.jl b/test/M_Boson_Fermion.jl
index 3cd5bc7d..7802b3a5 100644
--- a/test/M_Boson_Fermion.jl
+++ b/test/M_Boson_Fermion.jl
@@ -34,6 +34,8 @@
     @test nnz(L.data.A) == nnz(L(0)) == 43368
     L = addBosonDissipator(L, J)
     @test nnz(L.data.A) == nnz(L(0)) == 45590
+    @test isconstant(L)
+    @test iscached(L)
     ados = steadystate(L; verbose = false)
     @test ados.dims == L.dims
     @test length(ados) == L.N
diff --git a/test/M_Fermion.jl b/test/M_Fermion.jl
index dd8f4601..6cb4b174 100644
--- a/test/M_Fermion.jl
+++ b/test/M_Fermion.jl
@@ -33,6 +33,8 @@
     @test nnz(L.data.A) == nnz(L(0)) == 21318
     L = addFermionDissipator(L, J)
     @test nnz(L.data.A) == nnz(L(0)) == 22516
+    @test isconstant(L)
+    @test iscached(L)
     ados = steadystate(L; verbose = false)
     @test ados.dims == L.dims
     @test length(ados) == L.N
diff --git a/test/M_S.jl b/test/M_S.jl
index 10a54d65..a9c43308 100644
--- a/test/M_S.jl
+++ b/test/M_S.jl
@@ -12,6 +12,8 @@
     @test L.N == L_super.N == 1
     @test nnz(L.data.A) == nnz(L_super.data.A) == nnz(L(0)) == nnz(L_super(0)) == 8
     @test L.data == L_super.data
+    @test isconstant(L)
+    @test iscached(L)
     ados_list = heomsolve(L, ψ0, 0:1:t; reltol = 1e-8, abstol = 1e-10, progress_bar = Val(false)).ados
     ados = ados_list[end]
     @test ados.dims == L.dims
diff --git a/test/time_evolution.jl b/test/time_evolution.jl
index 21067c3c..80764156 100644
--- a/test/time_evolution.jl
+++ b/test/time_evolution.jl
@@ -41,10 +41,12 @@
     @test_throws ErrorException heomsolve(L, ados_wrong4, Δt, steps; progress_bar = Val(false))
 
     # using the method based on ODE solver
-    sol_e = heomsolve(L, ψ0, tlist; e_ops = e_ops, saveat = tlist, progress_bar = Val(true)) # also test progress bar
+    prob_e = HEOMsolveProblem(L, ψ0, tlist; e_ops = e_ops, saveat = tlist, progress_bar = Val(true)) # also test progress bar
+    sol_e = heomsolve(prob_e)
     sol_e2 = heomsolve(L, ψ0, tlist; e_ops = e_ops, progress_bar = Val(false))
     ρ_list_e = getRho.(sol_e.ados)
     expvals_e = sol_e.expect
+    @test !haskey(prob_e.prob.kwargs, :tstops) # tstops should not exist for time-independent cases
     @test show(devnull, MIME("text/plain"), sol_e) === nothing
     @test_logs (:warn,) evolution(L, ψ0, tlist; progress_bar = Val(false)) # deprecated function
     @test_throws ErrorException heomsolve(L, ψ0, tlist; verbose = true)
@@ -102,7 +104,7 @@
 
     tlist = 0:10:400
     p_fast = (amplitude = 0.5, delay = 20, integral = π / 2)
-    fastDD_sol = heomsolve(
+    fastDD_prob = HEOMsolveProblem(
         L,
         ψ0,
         tlist;
@@ -114,6 +116,7 @@
         abstol = 1e-12,
         progress_bar = Val(false),
     )
+    fastDD_sol = heomsolve(fastDD_prob)
     fastDD_ados = fastDD_sol.ados
     fastDD1 = real.(fastDD_sol.expect[1, :])
     fastDD2 = expect(P01, fastDD_ados)
@@ -161,6 +164,7 @@
         0.47479965067847246,
         0.47451220871416044,
     ]
+    @test fastDD_prob.prob.kwargs[:tstops] == tlist # tstops should be equal to tlist for time-dependent cases
     @test show(devnull, MIME("text/plain"), fastDD_sol) === nothing
     @test length(fastDD_sol.ados) == length(tlist)
     @test size(fastDD_sol.expect) == (1, length(tlist))