diff --git a/Project.toml b/Project.toml
index 067c9222..49f0ee0b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.14.1"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
+DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
@@ -22,6 +23,7 @@ SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
@@ -34,6 +36,7 @@ ArrayInterface = "6, 7"
 CUDA = "5"
 DiffEqBase = "6"
 DiffEqCallbacks = "2 - 3.1, 3.8"
+DiffEqNoiseProcess = "5"
 FFTW = "1.5"
 Graphs = "1.7"
 IncompleteLU = "0.2"
@@ -49,6 +52,7 @@ SciMLOperators = "0.3"
 SparseArrays = "<0.0.1, 1"
 SpecialFunctions = "2"
 StaticArraysCore = "1"
+StochasticDiffEq = "6"
 Test = "<0.0.1, 1"
 julia = "1.10"
 
diff --git a/docs/src/api.md b/docs/src/api.md
index 20c35bd2..bf0c796a 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -175,18 +175,22 @@ qeye
 ```@docs
 TimeEvolutionSol
 TimeEvolutionMCSol
+TimeEvolutionSSESol
 sesolveProblem
 mesolveProblem
-lr_mesolveProblem
 mcsolveProblem
 mcsolveEnsembleProblem
+ssesolveProblem
+ssesolveEnsembleProblem
+lr_mesolveProblem
 sesolve
 mesolve
-lr_mesolve
 mcsolve
+ssesolve
 dfd_mesolve
 dsf_mesolve
 dsf_mcsolve
+lr_mesolve
 liouvillian
 liouvillian_generalized
 steadystate
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index d6735644..1822ff09 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -22,18 +22,21 @@ import SciMLBase:
     remake,
     u_modified!,
     ODEProblem,
+    SDEProblem,
     EnsembleProblem,
     EnsembleThreads,
     FullSpecialize,
     CallbackSet,
     ContinuousCallback,
     DiscreteCallback
+import StochasticDiffEq: StochasticDiffEqAlgorithm, SRA1
 import SciMLOperators: MatrixOperator
 import LinearSolve: LinearProblem, SciMLLinearSolveAlgorithm, KrylovJL_MINRES, KrylovJL_GMRES
 import DiffEqBase: get_tstops
 import DiffEqCallbacks: PeriodicCallback, PresetTimeCallback, TerminateSteadyState
 import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm
 import OrdinaryDiffEqTsit5: Tsit5
+import DiffEqNoiseProcess: RealWienerProcess
 
 # other dependencies (in alphabetical order)
 import ArrayInterface: allowed_getindex, allowed_setindex!
@@ -73,6 +76,7 @@ include("time_evolution/mesolve.jl")
 include("time_evolution/lr_mesolve.jl")
 include("time_evolution/sesolve.jl")
 include("time_evolution/mcsolve.jl")
+include("time_evolution/ssesolve.jl")
 include("time_evolution/time_evolution_dynamical.jl")
 
 # Others
diff --git a/src/qobj/operator_sum.jl b/src/qobj/operator_sum.jl
index d19b4eac..909e27db 100644
--- a/src/qobj/operator_sum.jl
+++ b/src/qobj/operator_sum.jl
@@ -18,7 +18,7 @@ struct OperatorSum{CT<:Vector{<:Number},OT<:AbstractVector} <: AbstractQuantumOb
         mapreduce(x -> size(x) == size_1, &, operators) ||
             throw(DimensionMismatch("All the operators must have the same dimensions."))
         T = promote_type(
-            mapreduce(x -> eltype(x.data), promote_type, operators),
+            mapreduce(x -> eltype(x), promote_type, operators),
             mapreduce(eltype, promote_type, coefficients),
         )
         coefficients2 = T.(coefficients)
diff --git a/src/time_evolution/ssesolve.jl b/src/time_evolution/ssesolve.jl
new file mode 100644
index 00000000..d9fcbe34
--- /dev/null
+++ b/src/time_evolution/ssesolve.jl
@@ -0,0 +1,416 @@
+export ssesolveProblem, ssesolveEnsembleProblem, ssesolve
+
+#TODO: Check if works in GPU
+function _ssesolve_update_coefficients!(ψ, coefficients, sc_ops)
+    _get_en = op -> real(dot(ψ, op, ψ)) #this is en/2: <Sn + Sn'>/2 = Re<Sn>
+    @. coefficients[2:end-1] = _get_en(sc_ops) #coefficients of the OperatorSum: Σ Sn * en/2
+    coefficients[end] = -sum(x -> x^2, coefficients[2:end-1]) / 2 #this last coefficient is -Σen^2/8
+    return nothing
+end
+
+function ssesolve_drift!(du, u, p, t)
+    _ssesolve_update_coefficients!(u, p.K.coefficients, p.sc_ops)
+
+    mul!(du, p.K, u)
+
+    return nothing
+end
+
+function ssesolve_diffusion!(du, u, p, t)
+    @inbounds en = @view(p.K.coefficients[2:end-1])
+
+    # du:(H,W). du_reshaped:(H*W,). 
+    # H:Hilbert space dimension, W: number of sc_ops
+    du_reshaped = reshape(du, :)
+    mul!(du_reshaped, p.D, u) #du[:,i] = D[i] * u
+
+    du .-= u .* reshape(en, 1, :) #du[:,i] -= en[i] * u
+
+    return nothing
+end
+
+function _ssesolve_prob_func(prob, i, repeat)
+    internal_params = prob.p
+
+    noise = RealWienerProcess(
+        prob.tspan[1],
+        zeros(length(internal_params.sc_ops)),
+        zeros(length(internal_params.sc_ops)),
+        save_everystep = false,
+    )
+
+    noise_rate_prototype = similar(prob.u0, length(prob.u0), length(internal_params.sc_ops))
+
+    prm = merge(
+        internal_params,
+        (
+            expvals = similar(internal_params.expvals),
+            progr = ProgressBar(size(internal_params.expvals, 2), enable = false),
+        ),
+    )
+
+    return remake(prob, p = prm, noise = noise, noise_rate_prototype = noise_rate_prototype)
+end
+
+_ssesolve_output_func(sol, i) = (sol, false)
+
+function _ssesolve_generate_statistics!(sol, i, states, expvals_all)
+    sol_i = sol[:, i]
+    !isempty(sol_i.prob.kwargs[:saveat]) ?
+    states[i] = [QuantumObject(sol_i.u[i], dims = sol_i.prob.p.Hdims) for i in 1:length(sol_i.u)] : nothing
+
+    copyto!(view(expvals_all, i, :, :), sol_i.prob.p.expvals)
+    return nothing
+end
+
+@doc raw"""
+    ssesolveProblem(H::QuantumObject,
+        ψ0::QuantumObject,
+        tlist::AbstractVector;
+        sc_ops::Union{Nothing,AbstractVector}=nothing;
+        alg::StochasticDiffEqAlgorithm=SRA1()
+        e_ops::Union{Nothing,AbstractVector} = nothing,
+        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
+        params::NamedTuple=NamedTuple(),
+        progress_bar::Union{Val,Bool}=Val(true),
+        kwargs...)
+
+Generates the SDEProblem for the Stochastic Schrödinger time evolution of a quantum system. This is defined by the following stochastic differential equation:
+    
+    ```math
+    d|\psi(t)\rangle = -i K |\psi(t)\rangle dt + \sum_n M_n |\psi(t)\rangle dW_n(t)
+    ```
+
+where 
+    
+```math
+    K = \hat{H} + i \sum_n \left(\frac{e_j} C_n - \frac{1}{2} \sum_{j} C_n^\dagger C_n - \frac{e_j^2}{8}\right),
+    ```
+    ```math
+    M_n = C_n - \frac{e_n}{2},
+    ```
+    and
+    ```math
+    e_n = \langle C_n + C_n^\dagger \rangle.
+    ```
+
+Above, `C_n` is the `n`-th collapse operator and  `dW_j(t)` is the real Wiener increment associated to `C_n`.
+
+# Arguments
+
+- `H::QuantumObject`: The Hamiltonian of the system ``\hat{H}``.
+- `ψ0::QuantumObject`: The initial state of the system ``|\psi(0)\rangle``.
+- `tlist::AbstractVector`: The time list of the evolution.
+- `sc_ops::Union{Nothing,AbstractVector}=nothing`: List of stochastic collapse operators ``\{\hat{C}_n\}_n``.
+- `alg::StochasticDiffEqAlgorithm`: The algorithm used for the time evolution.
+- `e_ops::Union{Nothing,AbstractVector}=nothing`: The list of operators to be evaluated during the evolution.
+- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: The time-dependent Hamiltonian of the system. If `nothing`, the Hamiltonian is time-independent.
+- `params::NamedTuple`: The parameters of the system.
+- `progress_bar::Union{Val,Bool}`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities.
+- `kwargs...`: The keyword arguments passed to the `SDEProblem` constructor.
+
+# Notes
+
+- 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`). You can also specify `e_ops` and `saveat` separately.
+- The default tolerances in `kwargs` are given as `reltol=1e-2` and `abstol=1e-2`.
+- For more details about `alg` please refer to [`DifferentialEquations.jl` (SDE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/)
+- For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/)
+
+# Returns
+
+- `prob`: The `SDEProblem` for the Stochastic Schrödinger time evolution of the system.
+"""
+function ssesolveProblem(
+    H::QuantumObject{MT1,OperatorQuantumObject},
+    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+    tlist::AbstractVector,
+    sc_ops::Union{Nothing,AbstractVector} = nothing;
+    alg::StochasticDiffEqAlgorithm = SRA1(),
+    e_ops::Union{Nothing,AbstractVector} = nothing,
+    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
+    params::NamedTuple = NamedTuple(),
+    progress_bar::Union{Val,Bool} = Val(true),
+    kwargs...,
+) where {MT1<:AbstractMatrix,T2}
+    H.dims != ψ0.dims && throw(DimensionMismatch("The two quantum objects are not of the same Hilbert dimension."))
+
+    haskey(kwargs, :save_idxs) &&
+        throw(ArgumentError("The keyword argument \"save_idxs\" is not supported in QuantumToolbox."))
+
+    sc_ops isa Nothing &&
+        throw(ArgumentError("The list of collapse operators must be provided. Use sesolveProblem instead."))
+
+    !(H_t isa Nothing) && throw(ArgumentError("Time-dependent Hamiltonians are not currently supported in ssesolve."))
+
+    progress_bar_val = makeVal(progress_bar)
+
+    t_l = convert(Vector{Float64}, tlist) # Convert it into Float64 to avoid type instabilities for StochasticDiffEq.jl
+
+    ϕ0 = get_data(ψ0)
+
+    H_eff = get_data(H - T2(0.5im) * mapreduce(op -> op' * op, +, sc_ops))
+    sc_ops2 = get_data.(sc_ops)
+
+    coefficients = [1.0, fill(0.0, length(sc_ops) + 1)...]
+    operators = [-1im * H_eff, sc_ops2..., MT1(I(prod(H.dims)))]
+    K = OperatorSum(coefficients, operators)
+    _ssesolve_update_coefficients!(ϕ0, K.coefficients, sc_ops2)
+
+    D = reduce(vcat, sc_ops2)
+
+    progr = ProgressBar(length(t_l), enable = getVal(progress_bar_val))
+
+    if e_ops isa Nothing
+        expvals = Array{ComplexF64}(undef, 0, length(t_l))
+        e_ops2 = MT1[]
+        is_empty_e_ops = true
+    else
+        expvals = Array{ComplexF64}(undef, length(e_ops), length(t_l))
+        e_ops2 = get_data.(e_ops)
+        is_empty_e_ops = isempty(e_ops)
+    end
+
+    p = (
+        K = K,
+        D = D,
+        e_ops = e_ops2,
+        sc_ops = sc_ops2,
+        expvals = expvals,
+        progr = progr,
+        Hdims = H.dims,
+        H_t = H_t,
+        is_empty_e_ops = is_empty_e_ops,
+        params...,
+    )
+
+    saveat = e_ops isa Nothing ? t_l : [t_l[end]]
+    default_values = (DEFAULT_SDE_SOLVER_OPTIONS..., saveat = saveat)
+    kwargs2 = merge(default_values, kwargs)
+    kwargs3 = _generate_sesolve_kwargs(e_ops, progress_bar_val, t_l, kwargs2)
+
+    tspan = (t_l[1], t_l[end])
+    noise = RealWienerProcess(t_l[1], zeros(length(sc_ops)), zeros(length(sc_ops)), save_everystep = false)
+    noise_rate_prototype = similar(ϕ0, length(ϕ0), length(sc_ops))
+    return SDEProblem{true}(
+        ssesolve_drift!,
+        ssesolve_diffusion!,
+        ϕ0,
+        tspan,
+        p;
+        noise_rate_prototype = noise_rate_prototype,
+        noise = noise,
+        kwargs3...,
+    )
+end
+
+@doc raw"""
+    ssesolveEnsembleProblem(H::QuantumObject,
+        ψ0::QuantumObject,
+        tlist::AbstractVector;
+        sc_ops::Union{Nothing,AbstractVector} = nothing;
+        alg::StochasticDiffEqAlgorithm=SRA1()
+        e_ops::Union{Nothing,AbstractVector} = nothing,
+        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
+        params::NamedTuple=NamedTuple(),
+        prob_func::Function=_mcsolve_prob_func,
+        output_func::Function=_mcsolve_output_func,
+        kwargs...)
+
+Generates the SDE EnsembleProblem for the Stochastic Schrödinger time evolution of a quantum system. This is defined by the following stochastic differential equation:
+    
+    ```math
+    d|\psi(t)\rangle = -i K |\psi(t)\rangle dt + \sum_n M_n |\psi(t)\rangle dW_n(t)
+    ```
+
+where 
+    
+```math
+    K = \hat{H} + i \sum_n \left(\frac{e_j} C_n - \frac{1}{2} \sum_{j} C_n^\dagger C_n - \frac{e_j^2}{8}\right),
+    ```
+    ```math
+    M_n = C_n - \frac{e_n}{2},
+    ```
+    and
+    ```math
+    e_n = \langle C_n + C_n^\dagger \rangle.
+    ```
+
+Above, `C_n` is the `n`-th collapse operator and  `dW_j(t)` is the real Wiener increment associated to `C_n`.
+
+# Arguments
+
+- `H::QuantumObject`: The Hamiltonian of the system ``\hat{H}``.
+- `ψ0::QuantumObject`: The initial state of the system ``|\psi(0)\rangle``.
+- `tlist::AbstractVector`: The time list of the evolution.
+- `sc_ops::Union{Nothing,AbstractVector}=nothing`: List of stochastic collapse operators ``\{\hat{C}_n\}_n``.
+- `alg::StochasticDiffEqAlgorithm`: The algorithm used for the time evolution.
+- `e_ops::Union{Nothing,AbstractVector}=nothing`: The list of operators to be evaluated during the evolution.
+- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: The time-dependent Hamiltonian of the system. If `nothing`, the Hamiltonian is time-independent.
+- `params::NamedTuple`: The parameters of the system.
+- `prob_func::Function`: Function to use for generating the SDEProblem.
+- `output_func::Function`: Function to use for generating the output of a single trajectory.
+- `kwargs...`: The keyword arguments passed to the `SDEProblem` constructor.
+
+# Notes
+
+- 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`). You can also specify `e_ops` and `saveat` separately.
+- The default tolerances in `kwargs` are given as `reltol=1e-2` and `abstol=1e-2`.
+- For more details about `alg` please refer to [`DifferentialEquations.jl` (SDE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/)
+- For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/)
+
+# Returns
+
+- `prob::EnsembleProblem with SDEProblem`: The Ensemble SDEProblem for the Stochastic Shrödinger time evolution.
+"""
+function ssesolveEnsembleProblem(
+    H::QuantumObject{MT1,OperatorQuantumObject},
+    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+    tlist::AbstractVector,
+    sc_ops::Union{Nothing,AbstractVector} = nothing;
+    alg::StochasticDiffEqAlgorithm = SRA1(),
+    e_ops::Union{Nothing,AbstractVector} = nothing,
+    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
+    params::NamedTuple = NamedTuple(),
+    prob_func::Function = _ssesolve_prob_func,
+    output_func::Function = _ssesolve_output_func,
+    kwargs...,
+) where {MT1<:AbstractMatrix,T2}
+    prob_sse = ssesolveProblem(H, ψ0, tlist, sc_ops; alg = alg, e_ops = e_ops, H_t = H_t, params = params, kwargs...)
+
+    ensemble_prob = EnsembleProblem(prob_sse, prob_func = prob_func, output_func = output_func, safetycopy = false)
+
+    return ensemble_prob
+end
+
+@doc raw"""
+    ssesolve(H::QuantumObject,
+        ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+        tlist::AbstractVector,
+        sc_ops::Union{Nothing, AbstractVector}=nothing;
+        alg::StochasticDiffEqAlgorithm=SRA1(),
+        e_ops::Union{Nothing,AbstractVector}=nothing,
+        H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
+        params::NamedTuple=NamedTuple(),
+        n_traj::Int=1,
+        ensemble_method=EnsembleThreads(),
+        prob_func::Function=_mcsolve_prob_func,
+        output_func::Function=_mcsolve_output_func,
+        kwargs...)
+
+
+Stochastic Schrödinger equation evolution of a quantum system given the system Hamiltonian ``\hat{H}`` and a list of stochadtic collapse (jump) operators ``\{\hat{C}_n\}_n``.
+The stochastic evolution of the state ``|\psi(t)\rangle`` is defined by:
+    
+    ```math
+    d|\psi(t)\rangle = -i K |\psi(t)\rangle dt + \sum_n M_n |\psi(t)\rangle dW_n(t)
+    ```
+
+where 
+    
+```math
+    K = \hat{H} + i \sum_n \left(\frac{e_j} C_n - \frac{1}{2} \sum_{j} C_n^\dagger C_n - \frac{e_j^2}{8}\right),
+    ```
+    ```math
+    M_n = C_n - \frac{e_n}{2},
+    ```
+    and
+    ```math
+    e_n = \langle C_n + C_n^\dagger \rangle.
+    ```
+
+Above, `C_n` is the `n`-th collapse operator and  `dW_j(t)` is the real Wiener increment associated to `C_n`.
+
+
+# Arguments
+
+- `H::QuantumObject`: Hamiltonian of the system ``\hat{H}``.
+- `ψ0::QuantumObject`: Initial state of the system ``|\psi(0)\rangle``.
+- `tlist::AbstractVector`: List of times at which to save the state of the system.
+- `sc_ops::Union{Nothing,AbstractVector}=nothing`: List of stochastic collapse operators ``\{\hat{C}_n\}_n``.
+- `alg::StochasticDiffEqAlgorithm`: Algorithm to use for the time evolution.
+- `e_ops::Union{Nothing,AbstractVector}`: List of operators for which to calculate expectation values.
+- `H_t::Union{Nothing,Function,TimeDependentOperatorSum}`: Time-dependent part of the Hamiltonian.
+- `params::NamedTuple`: Dictionary of parameters to pass to the solver.
+- `seeds::Union{Nothing, Vector{Int}}`: List of seeds for the random number generator. Length must be equal to the number of trajectories provided.
+- `n_traj::Int`: Number of trajectories to use.
+- `ensemble_method`: Ensemble method to use.
+- `prob_func::Function`: Function to use for generating the SDEProblem.
+- `output_func::Function`: Function to use for generating the output of a single trajectory.
+- `kwargs...`: Additional keyword arguments to pass to the solver.
+
+# Notes
+
+- `ensemble_method` can be one of `EnsembleThreads()`, `EnsembleSerial()`, `EnsembleDistributed()`
+- 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`). You can also specify `e_ops` and `saveat` separately.
+- The default tolerances in `kwargs` are given as `reltol=1e-2` and `abstol=1e-2`.
+- For more details about `alg` please refer to [`DifferentialEquations.jl` (SDE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/)
+- For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/)
+
+# Returns
+
+- `sol::TimeEvolutionSSESol`: The solution of the time evolution. See also [`TimeEvolutionSSESol`](@ref)
+"""
+function ssesolve(
+    H::QuantumObject{MT1,OperatorQuantumObject},
+    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
+    tlist::AbstractVector,
+    sc_ops::Union{Nothing,AbstractVector} = nothing;
+    alg::StochasticDiffEqAlgorithm = SRA1(),
+    e_ops::Union{Nothing,AbstractVector} = nothing,
+    H_t::Union{Nothing,Function,TimeDependentOperatorSum} = nothing,
+    params::NamedTuple = NamedTuple(),
+    n_traj::Int = 1,
+    ensemble_method = EnsembleThreads(),
+    prob_func::Function = _ssesolve_prob_func,
+    output_func::Function = _ssesolve_output_func,
+    kwargs...,
+) where {MT1<:AbstractMatrix,T2}
+    ens_prob = ssesolveEnsembleProblem(
+        H,
+        ψ0,
+        tlist,
+        sc_ops;
+        alg = alg,
+        e_ops = e_ops,
+        H_t = H_t,
+        params = params,
+        prob_func = prob_func,
+        output_func = output_func,
+        kwargs...,
+    )
+
+    return ssesolve(ens_prob; alg = alg, n_traj = n_traj, ensemble_method = ensemble_method)
+end
+
+function ssesolve(
+    ens_prob::EnsembleProblem;
+    alg::StochasticDiffEqAlgorithm = SRA1(),
+    n_traj::Int = 1,
+    ensemble_method = EnsembleThreads(),
+)
+    sol = solve(ens_prob, alg, ensemble_method, trajectories = n_traj)
+    _sol_1 = sol[:, 1]
+
+    expvals_all = Array{ComplexF64}(undef, length(sol), size(_sol_1.prob.p.expvals)...)
+    states =
+        isempty(_sol_1.prob.kwargs[:saveat]) ? fill(QuantumObject[], length(sol)) :
+        Vector{Vector{QuantumObject}}(undef, length(sol))
+
+    foreach(i -> _ssesolve_generate_statistics!(sol, i, states, expvals_all), eachindex(sol))
+    expvals = dropdims(sum(expvals_all, dims = 1), dims = 1) ./ length(sol)
+
+    return TimeEvolutionSSESol(
+        n_traj,
+        _sol_1.t,
+        states,
+        expvals,
+        expvals_all,
+        sol.converged,
+        _sol_1.alg,
+        _sol_1.prob.kwargs[:abstol],
+        _sol_1.prob.kwargs[:reltol],
+    )
+end
diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl
index b0ba7671..bd60ede2 100644
--- a/src/time_evolution/time_evolution.jl
+++ b/src/time_evolution/time_evolution.jl
@@ -1,9 +1,10 @@
 export TimeDependentOperatorSum
-export TimeEvolutionSol, TimeEvolutionMCSol
+export TimeEvolutionSol, TimeEvolutionMCSol, TimeEvolutionSSESol
 
 export liouvillian, 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)
 
 @doc raw"""
     struct TimeEvolutionSol
@@ -95,6 +96,52 @@ function Base.show(io::IO, sol::TimeEvolutionMCSol)
     return nothing
 end
 
+@doc raw"""
+     struct TimeEvolutionSSESol
+ A structure storing the results and some information from solving trajectories of the Stochastic Shrodinger equation time evolution.
+ # Fields (Attributes)
+ - `n_traj::Int`: Number of trajectories
+ - `times::AbstractVector`: The time list of the evolution in each trajectory.
+ - `states::Vector{Vector{QuantumObject}}`: The list of result states in each trajectory.
+ - `expect::Matrix`: The expectation values (averaging all trajectories) corresponding to each time point in `times`.
+ - `expect_all::Array`: The expectation values corresponding to each trajectory and each time point in `times`
+ - `converged::Bool`: Whether the solution is converged or not.
+ - `alg`: The algorithm which is used during the solving process.
+ - `abstol::Real`: The absolute tolerance which is used during the solving process.
+ - `reltol::Real`: The relative tolerance which is used during the solving process.
+ """
+struct TimeEvolutionSSESol{
+    TT<:Vector{<:Real},
+    TS<:AbstractVector,
+    TE<:Matrix{ComplexF64},
+    TEA<:Array{ComplexF64,3},
+    T1<:Real,
+    T2<:Real,
+}
+    n_traj::Int
+    times::TT
+    states::TS
+    expect::TE
+    expect_all::TEA
+    converged::Bool
+    alg::StochasticDiffEqAlgorithm
+    abstol::T1
+    reltol::T2
+end
+
+function Base.show(io::IO, sol::TimeEvolutionSSESol)
+    print(io, "Solution of quantum trajectories\n")
+    print(io, "(converged: $(sol.converged))\n")
+    print(io, "--------------------------------\n")
+    print(io, "num_trajectories = $(sol.n_traj)\n")
+    print(io, "num_states = $(length(sol.states[1]))\n")
+    print(io, "num_expect = $(size(sol.expect, 1))\n")
+    print(io, "SDE alg.: $(sol.alg)\n")
+    print(io, "abstol = $(sol.abstol)\n")
+    print(io, "reltol = $(sol.reltol)\n")
+    return nothing
+end
+
 abstract type LindbladJumpCallbackType end
 
 struct ContinuousLindbladJumpCallback <: LindbladJumpCallbackType
diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl
index f86d0bfe..3d054698 100644
--- a/test/core-test/time_evolution.jl
+++ b/test/core-test/time_evolution.jl
@@ -42,7 +42,7 @@
         end
     end
 
-    @testset "mesolve and mcsolve" begin
+    @testset "mesolve, mcsolve, and ssesolve" begin
         N = 10
         a = destroy(N)
         a_d = a'
@@ -56,6 +56,7 @@
         sol_me3 = mesolve(H, psi0, t_l, c_ops, e_ops = e_ops, saveat = t_l, progress_bar = Val(false))
         sol_mc = mcsolve(H, psi0, t_l, c_ops, n_traj = 500, e_ops = e_ops, progress_bar = Val(false))
         sol_mc_states = mcsolve(H, psi0, t_l, c_ops, n_traj = 500, saveat = t_l, progress_bar = Val(false))
+        sol_sse = ssesolve(H, psi0, t_l, c_ops, n_traj = 500, e_ops = e_ops, progress_bar = Val(false))
 
         ρt_mc = [ket2dm.(normalize.(states)) for states in sol_mc_states.states]
         expect_mc_states = mapreduce(states -> expect.(Ref(e_ops[1]), states), hcat, ρt_mc)
@@ -63,8 +64,10 @@
 
         sol_me_string = sprint((t, s) -> show(t, "text/plain", s), sol_me)
         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 sum(abs.(sol_mc.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
+        @test sum(abs.(sol_sse.expect .- sol_me.expect)) / length(t_l) < 0.1
         @test length(sol_me.states) == 1
         @test size(sol_me.expect) == (length(e_ops), length(t_l))
         @test length(sol_me2.states) == length(t_l)
@@ -90,6 +93,16 @@
               "ODE alg.: $(sol_mc.alg)\n" *
               "abstol = $(sol_mc.abstol)\n" *
               "reltol = $(sol_mc.reltol)\n"
+        @test sol_sse_string ==
+              "Solution of quantum trajectories\n" *
+              "(converged: $(sol_sse.converged))\n" *
+              "--------------------------------\n" *
+              "num_trajectories = $(sol_sse.n_traj)\n" *
+              "num_states = $(length(sol_sse.states[1]))\n" *
+              "num_expect = $(size(sol_sse.expect, 1))\n" *
+              "SDE alg.: $(sol_sse.alg)\n" *
+              "abstol = $(sol_sse.abstol)\n" *
+              "reltol = $(sol_sse.reltol)\n"
 
         @testset "Type Inference mesolve" begin
             @inferred mesolveProblem(H, psi0, t_l, c_ops, e_ops = e_ops, progress_bar = Val(false))
@@ -115,6 +128,20 @@
             @inferred mcsolve(H, psi0, [0, 10], c_ops, n_traj = 500, progress_bar = Val(false))
             @inferred mcsolve(H, Qobj(zeros(Int64, N)), t_l, c_ops, n_traj = 500, progress_bar = Val(false))
         end
+
+        @testset "Type Inference ssesolve" begin
+            @inferred ssesolveEnsembleProblem(
+                H,
+                psi0,
+                t_l,
+                c_ops,
+                n_traj = 500,
+                e_ops = e_ops,
+                progress_bar = Val(false),
+            )
+            @inferred ssesolve(H, psi0, t_l, c_ops, n_traj = 500, e_ops = e_ops, progress_bar = Val(false))
+            @inferred ssesolve(H, psi0, t_l, c_ops, n_traj = 500, progress_bar = Val(true))
+        end
     end
 
     @testset "exceptions" begin