diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5df138717..d5c2c198f 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
 
 - Implement `EnrSpace` and corresponding functionality. ([#500])
 - Check for orthogonality breakdown in `Lanczos` solver for `spectrum`. ([#501])
+- Store both `times` and `times_states` in time evolution solutions. ([#506], [#504])
 - Fix errors in `Julia v1.12`. ([#507])
 
 ## [v0.32.1]
@@ -259,4 +260,6 @@ Release date: 2024-11-13
 [#494]: https://github.com/qutip/QuantumToolbox.jl/issues/494
 [#500]: https://github.com/qutip/QuantumToolbox.jl/issues/500
 [#501]: https://github.com/qutip/QuantumToolbox.jl/issues/501
+[#504]: https://github.com/qutip/QuantumToolbox.jl/issues/504
+[#506]: https://github.com/qutip/QuantumToolbox.jl/issues/506
 [#507]: https://github.com/qutip/QuantumToolbox.jl/issues/507
diff --git a/docs/src/users_guide/time_evolution/mcsolve.md b/docs/src/users_guide/time_evolution/mcsolve.md
index e5a3e0309..4f9aaa723 100644
--- a/docs/src/users_guide/time_evolution/mcsolve.md
+++ b/docs/src/users_guide/time_evolution/mcsolve.md
@@ -59,7 +59,7 @@ CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
 ```
 
 ```@example mcsolve
-times = LinRange(0.0, 10.0, 200)
+tlist = LinRange(0.0, 10.0, 200)
 
 ψ0 = tensor(fock(2, 0), fock(10, 8))
 a  = tensor(qeye(2), destroy(10))
@@ -69,7 +69,7 @@ H = 2 * π * a' * a + 2 * π * σm' * σm + 2 * π * 0.25 * (σm * a' + σm' * a
 c_ops = [sqrt(0.1) * a]
 e_ops = [a' * a, σm' * σm]
 
-sol_500 = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops)
+sol_500 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops)
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
@@ -78,8 +78,8 @@ ax = Axis(fig[1, 1],
     ylabel = "Expectation values",
     title = "Monte Carlo time evolution (500 trajectories)",
 )
-lines!(ax, times, real(sol_500.expect[1,:]), label = "cavity photon number", linestyle = :solid)
-lines!(ax, times, real(sol_500.expect[2,:]), label = "atom excitation probability", linestyle = :dash)
+lines!(ax, tlist, real(sol_500.expect[1,:]), label = "cavity photon number", linestyle = :solid)
+lines!(ax, tlist, real(sol_500.expect[2,:]), label = "atom excitation probability", linestyle = :dash)
 
 axislegend(ax, position = :rt)
 
@@ -93,9 +93,9 @@ We can also change the number of trajectories (`ntraj`). This can be used to exp
 ```@example mcsolve
 e_ops = [a' * a]
 
-sol_1   = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ntraj = 1)
-sol_10  = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ntraj = 10)
-sol_100 = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ntraj = 100)
+sol_1   = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj = 1)
+sol_10  = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj = 10)
+sol_100 = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj = 100)
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
@@ -104,9 +104,9 @@ ax = Axis(fig[1, 1],
     ylabel = "Expectation values",
     title = "Monte Carlo time evolution",
 )
-lines!(ax, times, real(sol_1.expect[1,:]), label = "1 trajectory", linestyle = :dashdot)
-lines!(ax, times, real(sol_10.expect[1,:]), label = "10 trajectories", linestyle = :dash)
-lines!(ax, times, real(sol_100.expect[1,:]), label = "100 trajectories", linestyle = :solid)
+lines!(ax, tlist, real(sol_1.expect[1,:]), label = "1 trajectory", linestyle = :dashdot)
+lines!(ax, tlist, real(sol_10.expect[1,:]), label = "10 trajectories", linestyle = :dash)
+lines!(ax, tlist, real(sol_100.expect[1,:]), label = "100 trajectories", linestyle = :solid)
 
 axislegend(ax, position = :rt)
 
@@ -126,8 +126,8 @@ Monte Carlo evolutions often need hundreds of trajectories to obtain sufficient
     See the [documentation of `DifferentialEquations.jl`](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/) for more details. Also, see Julia's documentation for more details about multithreading and adding more processes.
 
 ```julia
-sol_serial   = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ensemblealg=EnsembleSerial())
-sol_parallel = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ensemblealg=EnsembleThreads());
+sol_serial   = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ensemblealg=EnsembleSerial())
+sol_parallel = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ensemblealg=EnsembleThreads());
 ```
 
 !!! tip "Parallelization on a Cluster"
diff --git a/docs/src/users_guide/time_evolution/mesolve.md b/docs/src/users_guide/time_evolution/mesolve.md
index 9537d3bf8..dc5b3eb13 100644
--- a/docs/src/users_guide/time_evolution/mesolve.md
+++ b/docs/src/users_guide/time_evolution/mesolve.md
@@ -128,15 +128,14 @@ sol = mesolve(H, ψ0, tlist, c_ops, e_ops = [sigmaz(), sigmay()])
 We can therefore plot the expectation values:
 
 ```@example mesolve
-times = sol.times
 expt_z = real(sol.expect[1,:])
 expt_y = real(sol.expect[2,:])
 
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
-lines!(ax, times, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
-lines!(ax, times, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
+lines!(ax, tlist, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
+lines!(ax, tlist, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
 
 axislegend(ax, position = :rt)
 
@@ -210,8 +209,6 @@ tlist = LinRange(0.0, 10.0, 200)
 L = liouvillian(H_a + H_c + H_I, c_ops)
 sol = mesolve(L, ψ0, tlist, e_ops=[σm' * σm, a' * a])
 
-times = sol.times
-
 # expectation value of Number operator
 N_atom   = real(sol.expect[1,:])
 N_cavity = real(sol.expect[2,:])
@@ -219,8 +216,8 @@ N_cavity = real(sol.expect[2,:])
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
-lines!(ax, times, N_atom, label = "atom excitation probability", linestyle = :solid)
-lines!(ax, times, N_cavity, label = "cavity photon number", linestyle = :dash)
+lines!(ax, tlist, N_atom, label = "atom excitation probability", linestyle = :solid)
+lines!(ax, tlist, N_cavity, label = "cavity photon number", linestyle = :dash)
 
 axislegend(ax, position = :rt)
 
diff --git a/docs/src/users_guide/time_evolution/sesolve.md b/docs/src/users_guide/time_evolution/sesolve.md
index d3ba8a21c..70372a7f0 100644
--- a/docs/src/users_guide/time_evolution/sesolve.md
+++ b/docs/src/users_guide/time_evolution/sesolve.md
@@ -59,8 +59,8 @@ sol = sesolve(H, ψ0, tlist, e_ops = [sigmaz(), sigmay()])
     Here, we call [`sesolve`](@ref) directly instead of pre-defining [`sesolveProblem`](@ref) first (as shown previously).
 
 ```@example sesolve
-times = sol.times
-print(size(times))
+println(size(sol.times)) # time points corresponds to stored expectation values 
+println(size(sol.times_states)) # time points corresponds to stored states
 ```
 
 ```@example sesolve
@@ -77,8 +77,8 @@ expt_y = real(expt[2,:])
 # plot by CairoMakie.jl
 fig = Figure(size = (500, 350))
 ax = Axis(fig[1, 1], xlabel = "Time", ylabel = "Expectation values")
-lines!(ax, times, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
-lines!(ax, times, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
+lines!(ax, tlist, expt_z, label = L"\langle\hat{\sigma}_z\rangle", linestyle = :solid)
+lines!(ax, tlist, expt_y, label = L"\langle\hat{\sigma}_y\rangle", linestyle = :dash)
 
 axislegend(ax, position = :rb)
 
@@ -91,6 +91,9 @@ If the keyword argument `e_ops` is not specified (or given as an empty `Vector`)
 tlist = [0, 10]
 sol = sesolve(H, ψ0, tlist) # or specify: e_ops = []
 
+println(size(sol.times))
+println(size(sol.times_states))
+
 sol.states
 ```
 
@@ -104,9 +107,11 @@ sol = sesolve(H, ψ0, tlist, e_ops = [sigmay()], saveat = tlist)
 ```
 
 ```@example sesolve
+println(size(sol.times))
 sol.expect
 ```
 
 ```@example sesolve
+println(size(sol.times_states))
 sol.states
 ```
diff --git a/docs/src/users_guide/time_evolution/solution.md b/docs/src/users_guide/time_evolution/solution.md
index b6a694837..94609f18a 100644
--- a/docs/src/users_guide/time_evolution/solution.md
+++ b/docs/src/users_guide/time_evolution/solution.md
@@ -12,8 +12,9 @@ CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
 
 | **Fields (Attributes)** | **Description** |
 |:------------------------|:----------------|
-| `sol.times` | The time list of the evolution. |
-| `sol.states` | The list of result states. |
+| `sol.times` | The list of time points at which the expectation values are calculated during the evolution. |
+| `sol.times_states` | The list of time points at which the states are stored during the evolution. |
+| `sol.states` | The list of result states corresponding to each time point in `sol.times_states`. |
 | `sol.expect` | The expectation values corresponding to each time point in `sol.times`. |
 | `sol.alg` | The algorithm which is used during the solving process. |
 | `sol.abstol` | The absolute tolerance which is used during the solving process. |
@@ -54,7 +55,7 @@ nothing # hide
 
 Recall that `Julia` uses `Fortran`-style indexing that begins with one (i.e., `[1,:]` represents the 1-st observable, where `:` represents all values corresponding to `tlist`).
 
-Together with the array of times at which these expectation values are calculated:
+Together with the list of time points at which these expectation values are calculated:
 
 ```@example TE-solution
 times = sol.times
@@ -83,6 +84,13 @@ State vectors, or density matrices, are accessed in a similar manner:
 sol.states
 ```
 
+Together with the list of time points at which these states are stored:
+
+```@example TE-solution
+times = sol.times_states
+nothing # hide
+```
+
 Here, the solution contains only one (final) state. Because the `states` will be saved depend on the keyword argument `saveat` in `kwargs`. If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). One can also specify `e_ops` and `saveat` separately.
 
 Some other solvers can have other output.
diff --git a/src/correlations.jl b/src/correlations.jl
index 3eae092b5..a8ac92f27 100644
--- a/src/correlations.jl
+++ b/src/correlations.jl
@@ -19,7 +19,7 @@ end
         C::QuantumObject;
         kwargs...)
 
-Returns the two-times correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \hat{C}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
+Returns the two-time correlation function of three operators ``\hat{A}``, ``\hat{B}`` and ``\hat{C}``: ``\left\langle \hat{A}(t) \hat{B}(t + \tau) \hat{C}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
 
 If the initial state `ψ0` is given as `nothing`, then the [`steadystate`](@ref) will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
 """
@@ -96,7 +96,7 @@ end
         reverse::Bool=false,
         kwargs...)
 
-Returns the two-times correlation function of two operators ``\hat{A}`` and ``\hat{B}`` : ``\left\langle \hat{A}(t + \tau) \hat{B}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
+Returns the two-time correlation function of two operators ``\hat{A}`` and ``\hat{B}`` : ``\left\langle \hat{A}(t + \tau) \hat{B}(t) \right\rangle`` for a given initial state ``|\psi_0\rangle``.
 
 If the initial state `ψ0` is given as `nothing`, then the [`steadystate`](@ref) will be used as the initial state. Note that this is only implemented if `H` is constant ([`QuantumObject`](@ref)).
 
diff --git a/src/spectrum.jl b/src/spectrum.jl
index b97488575..9a6313ad4 100644
--- a/src/spectrum.jl
+++ b/src/spectrum.jl
@@ -313,7 +313,7 @@ end
 Calculate the power spectrum corresponding to a two-time correlation function using fast Fourier transform (FFT).
 
 # Parameters
-- `tlist::AbstractVector`: List of times at which the two-time correlation function is given.
+- `tlist::AbstractVector`: List of time points at which the two-time correlation function is given.
 - `corr::AbstractVector`: List of two-time correlations corresponding to the given time point in `tlist`.
 - `inverse::Bool`: Whether to use the inverse Fourier transform or not. Default to `false`.
 
diff --git a/src/time_evolution/brmesolve.jl b/src/time_evolution/brmesolve.jl
index cf2a15a25..85ffb2594 100644
--- a/src/time_evolution/brmesolve.jl
+++ b/src/time_evolution/brmesolve.jl
@@ -158,7 +158,7 @@ Solves for the dynamics of a system using the Bloch-Redfield master equation, gi
 
 - `H`: The system Hamiltonian. Must be an [`Operator`](@ref)
 - `ψ0`:  Initial state of the system $|\psi(0)\rangle$. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `a_ops`: Nested list with each element is a `Tuple` of operator-function pairs `(a_op, spectra)`, and the coupling [`Operator`](@ref) `a_op` must be hermitian with corresponding `spectra` being a `Function` of transition energy
 - `c_ops`: List of collapse operators corresponding to Lindblad dissipators
 - `sec_cutoff`: Cutoff for secular approximation. Use `-1` if secular approximation is not used when evaluating bath-coupling terms.
diff --git a/src/time_evolution/lr_mesolve.jl b/src/time_evolution/lr_mesolve.jl
index bc166a540..1c73010e1 100644
--- a/src/time_evolution/lr_mesolve.jl
+++ b/src/time_evolution/lr_mesolve.jl
@@ -7,10 +7,11 @@ A structure storing the results and some information from solving low-rank maste
 
 # Fields (Attributes)
 
-- `times::AbstractVector`: The time list of the evolution.
-- `states::Vector{QuantumObject}`: The list of result states.
+- `times::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
+- `times_states::AbstractVector`: The list of time points at which the states are stored during the evolution.
+- `states::Vector{QuantumObject}`: The list of result states corresponding to each time point in `times_states`.
 - `expect::Matrix`: The expectation values corresponding to each time point in `times`.
-- `fexpect::Matrix`: The function values at each time point.
+- `fexpect::Matrix`: The function values corresponding to each time point in `times`.
 - `retcode`: The return code from the solver.
 - `alg`: The algorithm which is used during the solving process.
 - `abstol::Real`: The absolute tolerance which is used during the solving process.
@@ -19,7 +20,8 @@ A structure storing the results and some information from solving low-rank maste
 - `B::Vector{QuantumObject}`: The `B` matrix of the low-rank algorithm at each time point.
 """
 struct TimeEvolutionLRSol{
-    TT<:AbstractVector{<:Real},
+    TT1<:AbstractVector{<:Real},
+    TT2<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Matrix{ComplexF64},
     RetT<:Enum,
@@ -29,7 +31,8 @@ struct TimeEvolutionLRSol{
     TSZB<:AbstractVector,
     TM<:Vector{<:Integer},
 }
-    times::TT
+    times::TT1
+    times_states::TT2
     states::TS
     expect::TE
     fexpect::TE
@@ -549,6 +552,7 @@ function lr_mesolve(prob::ODEProblem; kwargs...)
     ρt = map(x -> Qobj(x[1] * x[2] * x[1]', type = Operator(), dims = prob.p.Hdims), zip(zt, Bt))
 
     return TimeEvolutionLRSol(
+        prob.p.times,
         sol.t,
         ρt,
         prob.p.expvals,
diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl
index 7db7048ff..fe8da0876 100644
--- a/src/time_evolution/mcsolve.jl
+++ b/src/time_evolution/mcsolve.jl
@@ -89,7 +89,7 @@ If the environmental measurements register a quantum jump, the wave function und
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
@@ -195,7 +195,7 @@ If the environmental measurements register a quantum jump, the wave function und
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
@@ -320,7 +320,7 @@ If the environmental measurements register a quantum jump, the wave function und
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `alg`: The algorithm to use for the ODE solver. Default to `Tsit5()`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -412,6 +412,7 @@ function mcsolve(
     return TimeEvolutionMCSol(
         ntraj,
         ens_prob_mc.times,
+        _sol_1.t,
         states,
         expvals,
         expvals, # This is average_expect
diff --git a/src/time_evolution/mesolve.jl b/src/time_evolution/mesolve.jl
index 32305cb1b..85d6ed613 100644
--- a/src/time_evolution/mesolve.jl
+++ b/src/time_evolution/mesolve.jl
@@ -34,7 +34,7 @@ where
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
@@ -132,7 +132,7 @@ where
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `alg`: The algorithm for the ODE solver. The default value is `Tsit5()`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -206,6 +206,7 @@ function mesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit
 
     return TimeEvolutionSol(
         prob.times,
+        sol.t,
         ρt,
         _get_expvals(sol, SaveFuncMESolve),
         sol.retcode,
diff --git a/src/time_evolution/sesolve.jl b/src/time_evolution/sesolve.jl
index c5607f231..ef9d98f9b 100644
--- a/src/time_evolution/sesolve.jl
+++ b/src/time_evolution/sesolve.jl
@@ -30,7 +30,7 @@ Generate the ODEProblem for the Schrödinger time evolution of a quantum system:
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
 - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
@@ -104,7 +104,7 @@ Time evolution of a closed quantum system using the Schrödinger equation:
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `alg`: The algorithm for the ODE solver. The default is `Tsit5()`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used.
@@ -156,6 +156,7 @@ function sesolve(prob::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit
 
     return TimeEvolutionSol(
         prob.times,
+        sol.t,
         ψt,
         _get_expvals(sol, SaveFuncSESolve),
         sol.retcode,
diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl
index 710029197..aa6dbae9b 100644
--- a/src/time_evolution/smesolve.jl
+++ b/src/time_evolution/smesolve.jl
@@ -49,7 +49,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_i\}_i``. It can be either a `Vector` or a `Tuple`.
 - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -194,7 +194,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_i\}_i``. It can be either a `Vector` or a `Tuple`.
 - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -321,7 +321,7 @@ Above, ``\hat{C}_i`` represent the collapse operators related to pure dissipatio
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`Ket`](@ref), [`Operator`](@ref) or [`OperatorKet`](@ref).
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_i\}_i``. It can be either a `Vector` or a `Tuple`.
 - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided.
 - `alg`: The algorithm to use for the stochastic differential equation. Default is `SRIW1()` if `sc_ops` is an [`AbstractQuantumObject`](@ref) (diagonal noise), and `SRA2()` otherwise (non-diagonal noise).
@@ -425,6 +425,7 @@ function smesolve(
     return TimeEvolutionStochasticSol(
         ntraj,
         ens_prob.times,
+        _sol_1.t,
         states,
         expvals,
         expvals, # This is average_expect
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
index 6945dc81e..e68987c02 100644
--- a/src/time_evolution/ssesolve.jl
+++ b/src/time_evolution/ssesolve.jl
@@ -51,7 +51,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: `NullParameters` of parameters to pass to the solver.
@@ -189,7 +189,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and  ``dW_n(t)`` is t
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
 - `params`: `NullParameters` of parameters to pass to the solver.
@@ -317,7 +317,7 @@ Above, ``\hat{S}_n`` are the stochastic collapse operators and ``dW_n(t)`` is th
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
 - `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
-- `tlist`: List of times at which to save either the state or the expectation values of the system.
+- `tlist`: List of time points at which to save either the state or the expectation values of the system.
 - `sc_ops`: List of stochastic collapse operators ``\{\hat{S}_n\}_n``. It can be either a `Vector`, a `Tuple` or a [`AbstractQuantumObject`](@ref). It is recommended to use the last case when only one operator is provided.
 - `alg`: The algorithm to use for the stochastic differential equation. Default is `SRIW1()` if `sc_ops` is an [`AbstractQuantumObject`](@ref) (diagonal noise), and `SRA2()` otherwise (non-diagonal noise).
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -419,6 +419,7 @@ function ssesolve(
     return TimeEvolutionStochasticSol(
         ntraj,
         ens_prob.times,
+        _sol_1.t,
         states,
         expvals,
         expvals, # This is average_expect
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index 0116ff398..bfdd6942b 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -46,8 +46,9 @@ A structure storing the results and some information from solving time evolution
 
 # Fields (Attributes)
 
-- `times::AbstractVector`: The time list of the evolution.
-- `states::Vector{QuantumObject}`: The list of result states.
+- `times::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
+- `times_states::AbstractVector`: The list of time points at which the states are stored during the evolution.
+- `states::Vector{QuantumObject}`: The list of result states corresponding to each time point in `times_states`.
 - `expect::Union{AbstractMatrix,Nothing}`: The expectation values corresponding to each time point in `times`.
 - `retcode`: The return code from the solver.
 - `alg`: The algorithm which is used during the solving process.
@@ -55,7 +56,8 @@ A structure storing the results and some information from solving time evolution
 - `reltol::Real`: The relative tolerance which is used during the solving process.
 """
 struct TimeEvolutionSol{
-    TT<:AbstractVector{<:Real},
+    TT1<:AbstractVector{<:Real},
+    TT2<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Union{AbstractMatrix,Nothing},
     RETT<:Enum,
@@ -63,7 +65,8 @@ struct TimeEvolutionSol{
     AT<:Real,
     RT<:Real,
 }
-    times::TT
+    times::TT1
+    times_states::TT2
     states::TS
     expect::TE
     retcode::RETT
@@ -96,8 +99,9 @@ A structure storing the results and some information from solving quantum trajec
 # Fields (Attributes)
 
 - `ntraj::Int`: Number of trajectories
-- `times::AbstractVector`: The time list of the evolution.
-- `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory.
+- `times::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
+- `times_states::AbstractVector`: The list of time points at which the states are stored during the evolution.
+- `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times_states`.
 - `expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `times`.
 - `average_expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `times`.
 - `runs_expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `times`
@@ -109,7 +113,8 @@ A structure storing the results and some information from solving quantum trajec
 - `reltol::Real`: The relative tolerance which is used during the solving process.
 """
 struct TimeEvolutionMCSol{
-    TT<:AbstractVector{<:Real},
+    TT1<:AbstractVector{<:Real},
+    TT2<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Union{AbstractMatrix,Nothing},
     TEA<:Union{AbstractArray,Nothing},
@@ -120,7 +125,8 @@ struct TimeEvolutionMCSol{
     RT<:Real,
 }
     ntraj::Int
-    times::TT
+    times::TT1
+    times_states::TT2
     states::TS
     expect::TE
     average_expect::TE # Currently just a synonym for `expect`
@@ -158,8 +164,9 @@ A structure storing the results and some information from solving trajectories o
 # Fields (Attributes)
 
 - `ntraj::Int`: Number of trajectories
-- `times::AbstractVector`: The time list of the evolution.
-- `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory.
+- `times::AbstractVector`: The list of time points at which the expectation values are calculated during the evolution.
+- `times_states::AbstractVector`: The list of time points at which the states are stored during the evolution.
+- `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory and each time point in `times_states`.
 - `expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `times`.
 - `average_expect::Union{AbstractMatrix,Nothing}`: The expectation values (averaging all trajectories) corresponding to each time point in `times`.
 - `runs_expect::Union{AbstractArray,Nothing}`: The expectation values corresponding to each trajectory and each time point in `times`
@@ -169,7 +176,8 @@ A structure storing the results and some information from solving trajectories o
 - `reltol::Real`: The relative tolerance which is used during the solving process.
 """
 struct TimeEvolutionStochasticSol{
-    TT<:AbstractVector{<:Real},
+    TT1<:AbstractVector{<:Real},
+    TT2<:AbstractVector{<:Real},
     TS<:AbstractVector,
     TE<:Union{AbstractMatrix,Nothing},
     TEA<:Union{AbstractArray,Nothing},
@@ -179,7 +187,8 @@ struct TimeEvolutionStochasticSol{
     RT<:Real,
 }
     ntraj::Int
-    times::TT
+    times::TT1
+    times_states::TT2
     states::TS
     expect::TE
     average_expect::TE # Currently just a synonym for `expect`
diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl
index 2e07047da..98960764e 100644
--- a/src/time_evolution/time_evolution_dynamical.jl
+++ b/src/time_evolution/time_evolution_dynamical.jl
@@ -244,6 +244,7 @@ function dfd_mesolve(
 
     return TimeEvolutionSol(
         dfd_prob.times,
+        sol.t,
         ρt,
         _get_expvals(sol, SaveFuncMESolve),
         sol.retcode,
diff --git a/test/core-test/brmesolve.jl b/test/core-test/brmesolve.jl
index 68cfad8c1..63b6ba02c 100644
--- a/test/core-test/brmesolve.jl
+++ b/test/core-test/brmesolve.jl
@@ -104,11 +104,11 @@ end
     e_ops = pauli_vectors
     H = δ * 0.5 * sigmax() + ϵ * 0.5 * sigmaz()
     ψ0 = unit(2basis(2, 0) + basis(2, 1))
-    times = LinRange(0, 10, 100)
+    tlist = LinRange(0, 10, 100)
 
     for (me_c_ops, brme_c_ops, brme_a_ops) in arg_sets
-        me = mesolve(H, ψ0, times, me_c_ops, e_ops = e_ops, progress_bar = Val(false))
-        brme = brmesolve(H, ψ0, times, brme_a_ops, brme_c_ops, e_ops = e_ops, progress_bar = Val(false))
+        me = mesolve(H, ψ0, tlist, me_c_ops, e_ops = e_ops, progress_bar = Val(false))
+        brme = brmesolve(H, ψ0, tlist, brme_a_ops, brme_c_ops, e_ops = e_ops, progress_bar = Val(false))
 
         @test all(me.expect .== brme.expect)
     end
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index 883716405..e17db89c2 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -85,12 +85,15 @@ end
     @test prob.prob.f.f isa MatrixOperator
     @test sum(abs.(sol.expect[1, :] .- amp_rabi .* sin.(Ω_rabi * tlist) .^ 2)) / length(tlist) < 0.1
     @test length(sol.times) == length(tlist)
+    @test length(sol.times_states) == 1
     @test length(sol.states) == 1
     @test size(sol.expect) == (length(e_ops), length(tlist))
     @test length(sol2.times) == length(tlist)
+    @test length(sol2.times_states) == length(tlist)
     @test length(sol2.states) == length(tlist)
     @test sol2.expect === nothing
     @test length(sol3.times) == length(tlist)
+    @test length(sol3.times_states) == length(saveat)
     @test length(sol3.states) == length(saveat)
     @test size(sol3.expect) == (length(e_ops), length(tlist))
     @test sol.expect[1, saveat_idxs] ≈ expect(e_ops[1], sol3.states) atol = 1e-6
@@ -168,12 +171,15 @@ end
     @test TESetup.prob_me.prob.f.f isa MatrixOperator
     @test isket(sol_me5.states[1])
     @test length(sol_me.times) == length(tlist)
+    @test length(sol_me.times_states) == 1
     @test length(sol_me.states) == 1
     @test size(sol_me.expect) == (length(e_ops), length(tlist))
     @test length(sol_me2.times) == length(tlist)
+    @test length(sol_me2.times_states) == length(tlist)
     @test length(sol_me2.states) == length(tlist)
     @test sol_me2.expect === nothing
     @test length(sol_me3.times) == length(tlist)
+    @test length(sol_me3.times_states) == length(saveat)
     @test length(sol_me3.states) == length(saveat)
     @test size(sol_me3.expect) == (length(e_ops), length(tlist))
     @test sol_me3.expect[1, TESetup.saveat_idxs] ≈ expect(e_ops[1], sol_me3.states) atol = 1e-6
@@ -289,8 +295,10 @@ end
     @test sum(abs, vec(expect_mc_states_mean) .- vec(sol_me.expect[1, saveat_idxs])) / length(tlist) < 0.1
     @test sum(abs, vec(expect_mc_states_mean2) .- vec(sol_me.expect[1, saveat_idxs])) / length(tlist) < 0.1
     @test length(sol_mc.times) == length(tlist)
+    @test length(sol_mc.times_states) == 1
     @test size(sol_mc.expect) == (length(e_ops), length(tlist))
     @test length(sol_mc_states.times) == length(tlist)
+    @test length(sol_mc_states.times_states) == length(saveat)
     @test sol_mc_states.expect === nothing
 
     sol_mc_string = sprint((t, s) -> show(t, "text/plain", s), sol_mc)
@@ -391,6 +399,7 @@ end
 
     @test sum(abs, sol_sse.expect .- sol_me.expect) / length(tlist) < 0.1
     @test length(sol_sse.times) == length(tlist)
+    @test length(sol_sse.times_states) == 1
     @test size(sol_sse.expect) == (length(e_ops), length(tlist))
     @test isnothing(sol_sse.measurement)
     @test size(sol_sse2.measurement) == (length(c_ops), 20, length(tlist) - 1)
@@ -525,6 +534,7 @@ end
     @test sum(abs, sol_sme.expect .- sol_me.expect) / length(tlist) < 0.1
     @test sum(abs, sol_sme3.expect .- sol_me.expect) / length(tlist) < 0.1
     @test length(sol_sme.times) == length(tlist)
+    @test length(sol_sme.times_states) == 1
     @test size(sol_sme.expect) == (length(e_ops), length(tlist))
     @test isnothing(sol_sme.measurement)
     @test size(sol_sme2.measurement) == (length(sc_ops_sme), 20, length(tlist) - 1)