[Docs] update time evolution introduction and solution

docs/src/users_guide/time_evolution/intro.md

@@ -23,9 +23,12 @@ There are two kinds of quantum systems: open systems that interact with a larger
 
 The following table lists the solvers provided by `QuantumToolbox` for dynamic quantum systems and the corresponding type of solution returned by the solver:
 
-| **Equation** | **Function Call** | **Returned Solution** |
-|:-------------|:------------------|:----------------------|
-| Unitary evolution, Schrödinger equation | [`sesolve`](@ref) | [`TimeEvolutionSol`](@ref) |
-| Lindblad master eqn. or Von Neuman eqn. | [`mesolve`](@ref) | [`TimeEvolutionSol`](@ref) |
-| Monte Carlo evolution | [`mcsolve`](@ref) | [`TimeEvolutionMCSol`](@ref) |
-| Stochastic Schrödinger equation | [`ssesolve`](@ref) | [`TimeEvolutionSSESol`](@ref) |
+| **Equation** | **Function Call** | **Problem** | **Returned Solution** |
+|:-------------|:------------------|:------------|:----------------------|
+| Unitary evolution, Schrödinger equation | [`sesolve`](@ref) | [`sesolveProblem`](@ref) | [`TimeEvolutionSol`](@ref) |
+| Lindblad master eqn. or Von Neuman eqn.| [`mesolve`](@ref) | [`mesolveProblem`](@ref) | [`TimeEvolutionSol`](@ref) |
+| Monte Carlo evolution | [`mcsolve`](@ref) | [`mcsolveProblem`](@ref) [`mcsolveEnsembleProblem`](@ref) | [`TimeEvolutionMCSol`](@ref) |
+| Stochastic Schrödinger equation | [`ssesolve`](@ref) | [`ssesolveProblem`](@ref) [`ssesolveEnsembleProblem`](@ref) | [`TimeEvolutionSSESol`](@ref) |
+
+!!! note "Solving dynamics with pre-defined problems"
+    `QuantumToolbox` provides two different methods to solve the dynamics. One can use the function calls listed above by either taking all the operators (like Hamiltonian and collapse operators, etc.) as inputs directly, or generating the `prob`lems by yourself and take it as an input of the function call, e.g., `sesolve(prob)`.

docs/src/users_guide/time_evolution/solution.md

@@ -2,8 +2,86 @@
 
 ## [Time Evolution Solutions](@id doc-TE:Time-Evolution-Solutions)
 
-This page is still under construction, please visit [API](@ref doc-API) first.
+```@setup TE-solution
+using QuantumToolbox
+```
 
 ### Solution
+`QuantumToolbox` utilizes the powerful [`DifferentialEquation.jl`](https://docs.sciml.ai/DiffEqDocs/stable/) to simulate different kinds of quantum system dynamics. Thus, we will first look at the data structure used for returning the solution (`sol`) from [`DifferentialEquation.jl`](https://docs.sciml.ai/DiffEqDocs/stable/). The solution stores all the crucial data needed for analyzing and plotting the results of a simulation. A generic structure [`TimeEvolutionSol`](@ref) contains the following properties for storing simulation data:
 
-### Multiple trajectories solution
+| **Fields (Attributes)** | **Description** |
+|:------------------------|:----------------|
+| `sol.times` | The time list of the evolution. |
+| `sol.states` | The list of result 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. |
+| `sol.reltol` | The relative tolerance which is used during the solving process. |
+| `sol.retcode` (or `sol.converged`) | The returned status from the solver. |
+
+### Accessing data in solutions
+
+To understand how to access the data in solution, we will use an example as a guide, although we do not worry about the simulation details at this stage. The Schrödinger equation solver ([`sesolve`](@ref)) used in this example returns [`TimeEvolutionSol`](@ref):
+
+```@example TE-solution
+H = 0.5 * sigmax()
+ψ0 = basis(2, 0)
+e_ops = [
+    proj(basis(2, 0)),
+    proj(basis(2, 1)),
+    basis(2, 0) * basis(2, 1)'
+]
+tlist = LinRange(0, 10, 100)
+sol = sesolve(H, ψ0, tlist, e_ops = e_ops, progress_bar = Val(false)); nothing # hide
+```
+
+To see what is contained inside the solution, we can use the `print` function:
+
+```@example TE-solution
+print(sol)
+```
+
+It tells us the number of expectation values are computed and the number of states are stored. Now we have all the information needed to analyze the simulation results. To access the data for the three expectation values, one can do:
+
+```@example TE-solution
+expt1 = real(sol.expect[1,:])
+expt2 = real(sol.expect[2,:])
+expt3 = real(sol.expect[3,:]); 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:
+
+```@example TE-solution
+times = sol.times; nothing # hide
+```
+
+we can plot the resulting expectation values:
+
+```@example TE-solution
+using CairoMakie
+CairoMakie.enable_only_mime!(MIME"image/svg+xml"())
+
+fig = Figure()
+ax = Axis(fig[1, 1])
+lines!(ax, times, expt1, label = L"P_00")
+lines!(ax, times, expt2, label = L"P_11")
+lines!(ax, times, expt3, label = L"P_01")
+
+fig
+```
+
+State vectors, or density matrices, are accessed in a similar manner:
+
+```@example TE-solution
+sol.states
+```
+
+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 specified, the default value of `saveat=[tlist[end]]` (only save the final state), otherwise, `saveat=tlist` (saving the states corresponding to `tlist`). One can also specify `e_ops` and `saveat` separately.
+
+Some other solvers can have other output.
+
+### Multiple trajectories solution
+
+This part is still under construction, please visit [API](@ref doc-API) first.