diff --git a/CHANGELOG.md b/CHANGELOG.md
index 71fe85710..153eb036f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
 - Improve efficiency of `bloch_redfield_tensor` by avoiding unnecessary conversions. ([#509])
+- Support `SciMLOperators v1.4+`. ([#470])
 
 ## [v0.33.0]
 Release date: 2025-07-22
@@ -255,6 +256,7 @@ Release date: 2024-11-13
 [#455]: https://github.com/qutip/QuantumToolbox.jl/issues/455
 [#456]: https://github.com/qutip/QuantumToolbox.jl/issues/456
 [#460]: https://github.com/qutip/QuantumToolbox.jl/issues/460
+[#470]: https://github.com/qutip/QuantumToolbox.jl/issues/470
 [#472]: https://github.com/qutip/QuantumToolbox.jl/issues/472
 [#473]: https://github.com/qutip/QuantumToolbox.jl/issues/473
 [#476]: https://github.com/qutip/QuantumToolbox.jl/issues/476
diff --git a/Project.toml b/Project.toml
index 8b33fbb74..b93326f38 100644
--- a/Project.toml
+++ b/Project.toml
@@ -55,13 +55,13 @@ KernelAbstractions = "0.9.2"
 LaTeXStrings = "1.2"
 LinearAlgebra = "1"
 LinearSolve = "2, 3"
-Makie = "0.20, 0.21, 0.22, 0.23, 0.24"
+Makie = "0.20, 0.21, 0.22, 0.23"
 OrdinaryDiffEqCore = "1"
 OrdinaryDiffEqTsit5 = "1"
 Pkg = "1"
 Random = "1"
 SciMLBase = "2"
-SciMLOperators = "0.3, 0.4"
+SciMLOperators = "1.4"
 SparseArrays = "1"
 SpecialFunctions = "2"
 StaticArraysCore = "1"
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index a4d0b7841..afdbad0ed 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -16,6 +16,7 @@ import SciMLBase:
     u_modified!,
     NullParameters,
     ODEFunction,
+    SDEFunction,
     ODEProblem,
     SDEProblem,
     EnsembleProblem,
diff --git a/src/qobj/quantum_object_evo.jl b/src/qobj/quantum_object_evo.jl
index a4111f7c7..b7a165b4f 100644
--- a/src/qobj/quantum_object_evo.jl
+++ b/src/qobj/quantum_object_evo.jl
@@ -512,7 +512,7 @@ function (A::QuantumObjectEvolution)(
         )
     end
 
-    A.data(ψout.data, ψin.data, p, t)
+    A.data(ψout.data, ψin.data, nothing, p, t)
 
     return ψout
 end
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index bfdd6942b..aa68cfc00 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -395,7 +395,7 @@ end
 
 A struct to represent the diffusion operator. This is used to perform the diffusion process on N different Wiener processes.
 =#
-struct DiffusionOperator{T,OpType<:Tuple{Vararg{AbstractSciMLOperator}}} <: AbstractSciMLOperator{T}
+struct DiffusionOperator{T,OpType<:Tuple{Vararg{AbstractSciMLOperator}}}
     ops::OpType
     function DiffusionOperator(ops::OpType) where {OpType}
         T = mapreduce(eltype, promote_type, ops)
@@ -403,32 +403,29 @@ struct DiffusionOperator{T,OpType<:Tuple{Vararg{AbstractSciMLOperator}}} <: Abst
     end
 end
 
-@generated function update_coefficients!(L::DiffusionOperator, u, p, t)
-    ops_types = L.parameters[2].parameters
-    N = length(ops_types)
-    return quote
-        Base.@nexprs $N i -> begin
-            update_coefficients!(L.ops[i], u, p, t)
-        end
-
-        nothing
-    end
-end
-
-@generated function LinearAlgebra.mul!(v::AbstractVecOrMat, L::DiffusionOperator, u::AbstractVecOrMat)
+@generated function (L::DiffusionOperator)(w, v, p, t)
     ops_types = L.parameters[2].parameters
     N = length(ops_types)
     quote
-        M = length(u)
-        S = (size(v, 1), size(v, 2)) # This supports also `v` as a `Vector`
+        M = length(v)
+        S = (size(w, 1), size(w, 2)) # This supports also `w` as a `Vector`
         (S[1] == M && S[2] == $N) || throw(DimensionMismatch("The size of the output vector is incorrect."))
         Base.@nexprs $N i -> begin
-            mul!(@view(v[:, i]), L.ops[i], u)
+            op = L.ops[i]
+            op(@view(w[:, i]), v, v, p, t)
         end
-        return v
+        return w
     end
 end
 
+#TODO: Remove when a new release of SciMLBase.jl >2.104.0 is available
+(f::SDEFunction)(du, u, p, t) =
+    if f.f isa AbstractSciMLOperator
+        f.f(du, u, u, p, t)
+    else
+        f.f(du, u, p, t)
+    end
+
 #######################################
 
 function liouvillian_floquet(
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index e17db89c2..99bbb198c 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -484,7 +484,6 @@ end
     H = TESetup.H
     ψ0 = TESetup.ψ0
     tlist = TESetup.tlist
-    c_ops = TESetup.c_ops
     c_ops_sme = TESetup.c_ops_sme
     sc_ops_sme = TESetup.sc_ops_sme
     c_ops_sme2 = TESetup.c_ops_sme2