diff --git a/Project.toml b/Project.toml
index acf02c51..eb22b9b4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "HierarchicalEOM"
 uuid = "a62dbcb7-80f5-4d31-9a88-8b19fd92b128"
-version = "2.9.0"
+version = "2.10.0"
 authors = ["Yi-Te Huang"]
 
 [deps]
@@ -28,8 +28,8 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 HierarchicalEOM_CUDAExt = "CUDA"
 
 [compat]
-CUDA = "5"
-DiffEqCallbacks = "2 - 4"
+CUDA = "5.0 - 5.8, 5.9.4 - 5"
+DiffEqCallbacks = "4.2.1 - 4"
 FastExpm = "1.1"
 FillArrays = "1"
 IncompleteLU = "0.2"
@@ -42,7 +42,7 @@ ProgressMeter = "1.11.0"
 QuantumToolbox = "0.39"
 Reexport = "1"
 SciMLBase = "2.105"
-SciMLOperators = "1.4"
+SciMLOperators = "1.11"
 SparseArrays = "1"
 StaticArraysCore = "1"
 julia = "1.10"
diff --git a/src/steadystate.jl b/src/steadystate.jl
index 3c16c716..4408d485 100644
--- a/src/steadystate.jl
+++ b/src/steadystate.jl
@@ -4,12 +4,12 @@ Solve the steady state of the auxiliary density operators based on `LinearSolve.
 
 # Parameters
 - `M::AbstractHEOMLSMatrix` : the matrix given from HEOM model, where the parity should be `EVEN`.
-- `alg::SciMLLinearSolveAlgorithm` : The solving algorithm in package `LinearSolve.jl`. Default to `KrylovJL_GMRES(rtol=1e-12, atol=1e-14)`.
+- `alg::SciMLLinearSolveAlgorithm` : The solving algorithm in package `LinearSolve.jl`. Default to `KrylovJL_GMRES(rtol = 1e-12, atol = 1e-14)`.
 - `verbose::Bool` : To display verbose output or not. Defaults to `true`.
 - `kwargs` : The keyword arguments for the `LinearProblem`
 
 # Notes
-- For more details about `alg`, `kwargs`, and `LinearProblem`, please refer to [`LinearSolve.jl`](http://linearsolve.sciml.ai/stable/)
+- For more details about `alg`, `kwargs`, and `LinearProblem`, please refer to [`LinearSolve.jl`](http://linearsolve.sciml.ai/stable/).
 
 # Returns
 - `::ADOs` : The steady state of auxiliary density operators.
@@ -31,29 +31,14 @@ function QuantumToolbox.steadystate(
     A = _HandleSteadyStateMatrix(M)
     b = sparsevec([1], [1.0 + 0.0im], size(M, 1))
 
-    if verbose
-        println("Solving steady state for ADOs by linear-solve method...")
-        flush(stdout)
-    end
-    if (!haskey(kwargs, :Pl)) && (isa(A, SparseMatrixCSC))
-        if verbose
-            print("Calculating left preconditioner with ilu...")
-            flush(stdout)
-        end
-        kwargs = merge((; kwargs...), (Pl = ilu(A, τ = 0.01),))
-        if verbose
-            println("[DONE]")
-            flush(stdout)
-        end
-    end
-
     # solving x where A * x = b
     if verbose
-        print("Solving linear problem...")
+        print("Solving steady state for ADOs by linear-solve method...")
         flush(stdout)
     end
-    cache = init(LinearProblem(A, _HandleVectorType(M, b)), alg, kwargs...)
-    sol = solve!(cache)
+
+    prob = LinearProblem{true}(A, _HandleVectorType(M, b))
+    sol = solve(prob, alg; kwargs...)
     if verbose
         println("[DONE]")
         flush(stdout)