diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index 87d9535b..a7218d98 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -36,6 +36,10 @@ function _generate_sesolve_kwargs(e_ops, progress_bar::Val{false}, tlist, kwargs
     return _generate_sesolve_kwargs_with_callback(tlist, kwargs)
 end
 
+_sesolve_make_U_QobjEvo(H::QuantumObjectEvolution{<:MatrixOperator}) =
+    QobjEvo(MatrixOperator(-1im * H.data.A), dims = H.dims, type = Operator)
+_sesolve_make_U_QobjEvo(H) = QobjEvo(H, -1im)
+
 @doc raw"""
     sesolveProblem(
         H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
@@ -88,7 +92,7 @@ function sesolveProblem(
 
     tlist = convert(Vector{_FType(ψ0)}, tlist) # Convert it to support GPUs and avoid type instabilities for OrdinaryDiffEq.jl
 
-    H_evo = QobjEvo(H, -1im) # pre-multiply by -i
+    H_evo = _sesolve_make_U_QobjEvo(H) # Multiply by -i
     isoper(H_evo) || throw(ArgumentError("The Hamiltonian must be an Operator."))
     check_dims(H_evo, ψ0)
 
diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl
index 8830de40..aa9660b0 100644
--- a/src/time_evolution/time_evolution_dynamical.jl
+++ b/src/time_evolution/time_evolution_dynamical.jl
@@ -564,8 +564,8 @@ function _DSF_mcsolve_Affect!(integrator)
     @. e_ops0 = get_data(e_ops2)
     @. c_ops0 = get_data(c_ops2)
     H_nh = lmul!(convert(eltype(ψt), 0.5im), mapreduce(op -> op' * op, +, c_ops0))
-    # By doing this, we are assuming that the system is time-independent and f is a ScaledOperator
-    copyto!(integrator.f.f.L.A, H(op_l2, dsf_params).data - H_nh)
+    # By doing this, we are assuming that the system is time-independent and f is a MatrixOperator
+    copyto!(integrator.f.f.A, lmul!(-1im, H(op_l2, dsf_params).data - H_nh))
     return u_modified!(integrator, true)
 end
 
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 64c42a6a..83adb961 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -98,7 +98,7 @@
         sol_mc_string = sprint((t, s) -> show(t, "text/plain", s), sol_mc)
         sol_sse_string = sprint((t, s) -> show(t, "text/plain", s), sol_sse)
         @test prob_me.f.f isa MatrixOperator
-        @test prob_mc.f.f isa ScaledOperator # TODO: can be optimized as MatrixOperator
+        @test prob_mc.f.f isa MatrixOperator
         @test sum(abs.(sol_mc.expect .- sol_me.expect)) / length(t_l) < 0.1
         @test sum(abs.(sol_mc2.expect .- sol_me.expect)) / length(t_l) < 0.1
         @test sum(abs.(vec(expect_mc_states_mean) .- vec(sol_me.expect))) / length(t_l) < 0.1