diff --git a/CHANGELOG.md b/CHANGELOG.md
index 02bbec4fe..f829b30f7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Change default solver detection in `eigensolve` when using `sigma` keyword argument (shift-inverse algorithm). If the operator is a `SparseMatrixCSC`, the default solver is `UMFPACKFactorization`, otherwise it is automatically chosen by LinearSolve.jl, depending on the type of the operator. ([#580]) 
 - Add keyword argument `assume_hermitian` to `liouvillian`. This allows users to disable the assumption that the Hamiltonian is Hermitian. ([#581])
+- Use LinearSolve's internal methods for preconditioners in `SteadyStateLinearSolver`. ([#588])
 
 ## [v0.38.1]
 Release date: 2025-10-27
@@ -362,3 +363,4 @@ Release date: 2024-11-13
 [#579]: https://github.com/qutip/QuantumToolbox.jl/issues/579
 [#580]: https://github.com/qutip/QuantumToolbox.jl/issues/580
 [#581]: https://github.com/qutip/QuantumToolbox.jl/issues/581
+[#588]: https://github.com/qutip/QuantumToolbox.jl/issues/588
diff --git a/docs/src/users_guide/steadystate.md b/docs/src/users_guide/steadystate.md
index a74353b90..fac4f4083 100644
--- a/docs/src/users_guide/steadystate.md
+++ b/docs/src/users_guide/steadystate.md
@@ -39,10 +39,11 @@ To change a solver, use the keyword argument `solver`, for example:
 
 Although it is not obvious, the [`SteadyStateDirectSolver()`](@ref SteadyStateDirectSolver) and [`SteadyStateEigenSolver()`](@ref SteadyStateEigenSolver) methods all use an LU decomposition internally and thus can have a large memory overhead. In contrast, for [`SteadyStateLinearSolver()`](@ref SteadyStateLinearSolver), iterative algorithms provided by [`LinearSolve.jl`](https://docs.sciml.ai/LinearSolve/stable/solvers/solvers/), such as `KrylovJL_GMRES()` and `KrylovJL_BICGSTAB()`, do not factor the matrix and thus take less memory than the LU methods and allow, in principle, for extremely large system sizes. The downside is that these methods can take much longer than the direct method as the condition number of the Liouvillian matrix is large, indicating that these iterative methods require a large number of iterations for convergence. 
 
-To overcome this, one can provide preconditioner that solves for an approximate inverse for the (modified) Liouvillian, thus better conditioning the problem, leading to faster convergence. The left and right preconditioner can be specified as the keyword argument `Pl` and `Pr`, respectively:
+To overcome this, preconditioners can be defined directly in the solver algorithm using LinearSolve's precs parameter. For example:
 ```julia
-solver = SteadyStateLinearSolver(alg = KrylovJL_GMRES(), Pl = Pl, Pr = Pr)
+SteadyStateLinearSolver(alg = KrylovJL_GMRES(; precs = (A, p) -> (ilu(A, τ = 0.01), I)))
 ```
+
 The use of a preconditioner can actually make these iterative methods faster than the other solution methods. The problem with precondioning is that it is only well defined for Hermitian matrices. Since the Liouvillian is non-Hermitian, the ability to find a good preconditioner is not guaranteed. Moreover, if a preconditioner is found, it is not guaranteed to have a good condition number. 
 
 Furthermore, `QuantiumToolbox` can take advantage of the Intel (Pardiso) LU solver in the Intel Math Kernel library (Intel-MKL) by using [`LinearSolve.jl`](https://docs.sciml.ai/LinearSolve/stable/) and either [`Pardiso.jl`](https://github.com/JuliaSparse/Pardiso.jl) or [`MKL_jll.jl`](https://github.com/JuliaBinaryWrappers/MKL_jll.jl):
diff --git a/src/steadystate.jl b/src/steadystate.jl
index 81ccd6833..edbe28c97 100644
--- a/src/steadystate.jl
+++ b/src/steadystate.jl
@@ -23,23 +23,20 @@ A solver which solves [`steadystate`](@ref) by finding the zero (or lowest) eige
 struct SteadyStateEigenSolver <: SteadyStateSolver end
 
 @doc raw"""
-    SteadyStateLinearSolver(alg = KrylovJL_GMRES(), Pl = nothing, Pr = nothing)
+    SteadyStateLinearSolver(
+        alg = KrylovJL_GMRES(; precs = (A, p) -> A isa SparseMatrixCSC ? (ilu(A, τ = 0.01), I) : (I, I))
+    )
 
 A solver which solves [`steadystate`](@ref) by finding the inverse of Liouvillian [`SuperOperator`](@ref) using the `alg`orithms given in [`LinearSolve.jl`](https://docs.sciml.ai/LinearSolve/stable/).
 
 # Arguments
 - `alg::SciMLLinearSolveAlgorithm=KrylovJL_GMRES()`: algorithms given in [`LinearSolve.jl`](https://docs.sciml.ai/LinearSolve/stable/)
-- `Pl::Union{Function,Nothing}=nothing`: left preconditioner, see documentation [Solving for Steady-State Solutions](@ref doc:Solving-for-Steady-State-Solutions) for more details.
-- `Pr::Union{Function,Nothing}=nothing`: right preconditioner, see documentation [Solving for Steady-State Solutions](@ref doc:Solving-for-Steady-State-Solutions) for more details.
+
+# Note
+Refer to [`LinearSolve.jl`](https://docs.sciml.ai/LinearSolve/stable/) for more details about the available algorithms. For example, the preconditioners can be defined directly in the solver like: `SteadyStateLinearSolver(alg = KrylovJL_GMRES(; precs = (A, p) -> (I, Diagonal(A))))`.
 """
-Base.@kwdef struct SteadyStateLinearSolver{
-    MT<:Union{SciMLLinearSolveAlgorithm,Nothing},
-    PlT<:Union{Function,Nothing},
-    PrT<:Union{Function,Nothing},
-} <: SteadyStateSolver
-    alg::MT = KrylovJL_GMRES()
-    Pl::PlT = nothing
-    Pr::PrT = nothing
+Base.@kwdef struct SteadyStateLinearSolver{MT<:Union{SciMLLinearSolveAlgorithm,Nothing}} <: SteadyStateSolver
+    alg::MT = KrylovJL_GMRES(; precs = (A, p) -> A isa SparseMatrixCSC ? (ilu(A, τ = 0.01), I) : (I, I))
 end
 
 @doc raw"""
@@ -159,14 +156,8 @@ function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateLinear
     L_tmp = L_tmp + Tn
 
     (haskey(kwargs, :Pl) || haskey(kwargs, :Pr)) && error("The use of preconditioners must be defined in the solver.")
-    if !isnothing(solver.Pl)
-        kwargs = merge((; kwargs...), (Pl = solver.Pl(L_tmp),))
-    elseif isa(L_tmp, SparseMatrixCSC)
-        kwargs = merge((; kwargs...), (Pl = ilu(L_tmp, τ = 0.01),))
-    end
-    !isnothing(solver.Pr) && (kwargs = merge((; kwargs...), (Pr = solver.Pr(L_tmp),)))
 
-    prob = LinearProblem(L_tmp, v0)
+    prob = LinearProblem{true}(L_tmp, v0)
     ρss_vec = solve(prob, solver.alg; kwargs...).u
 
     ρss = reshape(ρss_vec, N, N)
@@ -407,16 +398,10 @@ function _steadystate_fourier(
     v0[n_max*N+1] = weight
 
     (haskey(kwargs, :Pl) || haskey(kwargs, :Pr)) && error("The use of preconditioners must be defined in the solver.")
-    if !isnothing(solver.Pl)
-        kwargs = merge((; kwargs...), (Pl = solver.Pl(M),))
-    elseif isa(M, SparseMatrixCSC)
-        kwargs = merge((; kwargs...), (Pl = ilu(M, τ = 0.01),))
-    end
-    !isnothing(solver.Pr) && (kwargs = merge((; kwargs...), (Pr = solver.Pr(M),)))
     !haskey(kwargs, :abstol) && (kwargs = merge((; kwargs...), (abstol = tol,)))
     !haskey(kwargs, :reltol) && (kwargs = merge((; kwargs...), (reltol = tol,)))
 
-    prob = LinearProblem(M, v0)
+    prob = LinearProblem{true}(M, v0)
     ρtot = solve(prob, solver.alg; kwargs...).u
 
     offset1 = n_max * N