diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2439ab0a6..2c24b137e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Change save callbacks from `PresetTimeCallback` to `FunctionCallingCallback`. ([#410])
 - Align `eigenstates` and `eigenenergies` to QuTiP. ([#411])
 - Introduce `vector_to_operator` and `operator_to_vector`. ([#413])
+- Change structure of tolerances between pure and mixed state dynamics. ([#415])
 
 ## [v0.27.0]
 Release date: 2025-02-14
@@ -149,3 +150,4 @@ Release date: 2024-11-13
 [#410]: https://github.com/qutip/QuantumToolbox.jl/issues/410
 [#411]: https://github.com/qutip/QuantumToolbox.jl/issues/411
 [#413]: https://github.com/qutip/QuantumToolbox.jl/issues/413
+[#415]: https://github.com/qutip/QuantumToolbox.jl/issues/415
diff --git a/src/time_evolution/lr_mesolve.jl b/src/time_evolution/lr_mesolve.jl
index 873e7e9ed..cb35f148d 100644
--- a/src/time_evolution/lr_mesolve.jl
+++ b/src/time_evolution/lr_mesolve.jl
@@ -450,7 +450,7 @@ function lr_mesolveProblem(
     p.Si .= pinv(Hermitian(p.S), atol = opt.atol_inv)
 
     # Initialization of ODEProblem's kwargs
-    default_values = (DEFAULT_ODE_SOLVER_OPTIONS..., saveat = [t_l[end]])
+    default_values = (DEFAULT_ODE_ME_SOLVER_OPTIONS..., saveat = [t_l[end]])
     kwargs2 = merge(default_values, kwargs)
 
     # Initialization of Callbacks
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index 210cf37a7..4856dd6e9 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -132,7 +132,7 @@ function mcsolveProblem(
     T = Base.promote_eltype(H_eff_evo, ψ0)
 
     # We disable the progress bar of the sesolveProblem because we use a global progress bar for all the trajectories
-    default_values = (DEFAULT_ODE_SOLVER_OPTIONS..., progress_bar = Val(false))
+    default_values = (DEFAULT_ODE_SE_MC_SOLVER_OPTIONS..., progress_bar = Val(false))
     kwargs2 = _merge_saveat(tlist, e_ops, default_values; kwargs...)
     kwargs3 = _generate_mcsolve_kwargs(ψ0, T, e_ops, tlist, c_ops, jump_callback, rng, kwargs2)
 
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index ee78e7eac..70da8afc5 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -79,7 +79,7 @@ function mesolveProblem(
     ρ0 = to_dense(_CType(T), mat2vec(ket2dm(ψ0).data)) # Convert it to dense vector with complex element type
     L = L_evo.data
 
-    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...)
+    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_ME_SOLVER_OPTIONS; kwargs...)
     kwargs3 = _generate_se_me_kwargs(e_ops, makeVal(progress_bar), tlist, kwargs2, SaveFuncMESolve)
 
     tspan = (tlist[1], tlist[end])
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index 571b1a396..6f611eb4e 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -69,7 +69,7 @@ function sesolveProblem(
     ψ0 = to_dense(_CType(T), get_data(ψ0)) # Convert it to dense vector with complex element type
     U = H_evo.data
 
-    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...)
+    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_SE_MC_SOLVER_OPTIONS; kwargs...)
     kwargs3 = _generate_se_me_kwargs(e_ops, makeVal(progress_bar), tlist, kwargs2, SaveFuncSESolve)
 
     tspan = (tlist[1], tlist[end])
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
index 4890030d2..c36c11634 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -118,7 +118,7 @@ function smesolveProblem(
     end
     D = DiffusionOperator(D_l)
 
-    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_SDE_SOLVER_OPTIONS; kwargs...)
+    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_SDE_ME_SOLVER_OPTIONS; kwargs...)
     kwargs3 = _generate_stochastic_kwargs(
         e_ops,
         sc_ops_list,
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index 5d7e89d0b..7d82a52f9 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -118,7 +118,7 @@ function ssesolveProblem(
     D_l = map(op -> op + _ScalarOperator_e(op, -) * IdentityOperator(prod(dims)), sc_ops_evo_data)
     D = DiffusionOperator(D_l)
 
-    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_SDE_SOLVER_OPTIONS; kwargs...)
+    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_SDE_SE_SOLVER_OPTIONS; kwargs...)
     kwargs3 = _generate_stochastic_kwargs(
         e_ops,
         sc_ops_list,
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index b487d6511..fc4534825 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -2,8 +2,11 @@ export TimeEvolutionSol, TimeEvolutionMCSol, TimeEvolutionStochasticSol
 
 export liouvillian_floquet, liouvillian_generalized
 
-const DEFAULT_ODE_SOLVER_OPTIONS = (abstol = 1e-8, reltol = 1e-6, save_everystep = false, save_end = true)
-const DEFAULT_SDE_SOLVER_OPTIONS = (abstol = 1e-2, reltol = 1e-2, save_everystep = false, save_end = true)
+# Dynamics of density matrices generally support larger tolerances
+const DEFAULT_ODE_SE_MC_SOLVER_OPTIONS = (abstol = 1e-8, reltol = 1e-6, save_everystep = false, save_end = true)
+const DEFAULT_ODE_ME_SOLVER_OPTIONS = (abstol = 1e-7, reltol = 1e-5, save_everystep = false, save_end = true)
+const DEFAULT_SDE_SE_SOLVER_OPTIONS = (abstol = 5e-3, reltol = 1e-2, save_everystep = false, save_end = true)
+const DEFAULT_SDE_ME_SOLVER_OPTIONS = (abstol = 1e-2, reltol = 1e-2, save_everystep = false, save_end = true)
 const COL_TIMES_WHICH_INIT_SIZE = 200
 
 @doc raw"""