diff --git a/CHANGELOG.md b/CHANGELOG.md
index a1e7269f1..f71b63114 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
+- Add support of `QobjEvo` for `steadystate` (ODE solver only). ([#536])
+
 ## [v0.34.1]
 Release date: 2025-08-23
 
@@ -301,3 +303,4 @@ Release date: 2024-11-13
 [#517]: https://github.com/qutip/QuantumToolbox.jl/issues/517
 [#520]: https://github.com/qutip/QuantumToolbox.jl/issues/520
 [#531]: https://github.com/qutip/QuantumToolbox.jl/issues/531
+[#536]: https://github.com/qutip/QuantumToolbox.jl/issues/536
diff --git a/src/steadystate.jl b/src/steadystate.jl
index a81fe3b3f..7e55911ac 100644
--- a/src/steadystate.jl
+++ b/src/steadystate.jl
@@ -96,7 +96,7 @@ end
 
 @doc raw"""
     steadystate(
-        H::QuantumObject{OpType},
+        H::AbstractQuantumObject{OpType},
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         solver::SteadyStateSolver = SteadyStateDirectSolver(),
         kwargs...,
@@ -111,7 +111,7 @@ Solve the stationary state based on different solvers.
 - `kwargs`: The keyword arguments for the solver.
 """
 function steadystate(
-    H::QuantumObject{OpType},
+    H::AbstractQuantumObject{OpType},
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     solver::SteadyStateSolver = SteadyStateDirectSolver(),
     kwargs...,
@@ -199,7 +199,7 @@ function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateDirect
     return QuantumObject(ρss, Operator(), L.dimensions)
 end
 
-function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateODESolver; kwargs...)
+function _steadystate(L::AbstractQuantumObject{SuperOperator}, solver::SteadyStateODESolver; kwargs...)
     tmax = solver.tmax
 
     ψ0 = isnothing(solver.ψ0) ? rand_ket(L.dimensions) : solver.ψ0
@@ -212,17 +212,26 @@ function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateODESol
     sol = mesolve(
         L,
         ψ0,
-        [ftype(0), ftype(tmax)],
+        [ftype(0), ftype(tmax)];
+        alg = solver.alg,
         progress_bar = Val(false),
         save_everystep = false,
         saveat = ftype[],
         callback = cb,
+        kwargs...,
     )
 
     ρss = sol.states[end]
     return ρss
 end
 
+_steadystate(
+    L::QuantumObjectEvolution{SuperOperator},
+    solver::T;
+    kwargs...,
+) where {T<:Union{SteadyStateDirectSolver,SteadyStateEigenSolver,SteadyStateLinearSolver}} =
+    throw(ArgumentError("$(get_typename_wrapper(solver)) does not support QobjEvo."))
+
 struct SteadyStateODECondition{CT<:AbstractArray}
     cache::CT
 end
diff --git a/test/core-test/steady_state.jl b/test/core-test/steady_state.jl
index d78c93c33..6151af349 100644
--- a/test/core-test/steady_state.jl
+++ b/test/core-test/steady_state.jl
@@ -3,6 +3,7 @@
     a = destroy(N)
     a_d = a'
     H = a_d * a + 0.1 * (a + a_d)
+    Ht = QobjEvo(H, (p, t) -> 1) # for test throw
     c_ops = [sqrt(0.1) * a]
     e_ops = [a_d * a]
     psi0 = fock(N, 3)
@@ -17,18 +18,17 @@
     solver = SteadyStateDirectSolver()
     ρ_ss = steadystate(H, c_ops, solver = solver)
     @test tracedist(rho_me, ρ_ss) < 1e-4
+    @test_throws ArgumentError steadystate(Ht, c_ops, solver = solver)
 
     solver = SteadyStateLinearSolver()
     ρ_ss = steadystate(H, c_ops, solver = solver)
     @test tracedist(rho_me, ρ_ss) < 1e-4
-
-    solver = SteadyStateLinearSolver()
-    ρ_ss = steadystate(H, c_ops, solver = solver)
-    @test tracedist(rho_me, ρ_ss) < 1e-4
+    @test_throws ArgumentError steadystate(Ht, c_ops, solver = solver)
 
     solver = SteadyStateEigenSolver()
     ρ_ss = steadystate(H, c_ops, solver = solver)
     @test tracedist(rho_me, ρ_ss) < 1e-4
+    @test_throws ArgumentError steadystate(Ht, c_ops, solver = solver)
 
     @testset "Type Inference (steadystate)" begin
         L = liouvillian(H, c_ops)