diff --git a/.gitignore b/.gitignore
index b07886a94..014c79889 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,4 +10,8 @@ Manifest.toml
 benchmarks/benchmarks_output.json
 
 .ipynb_checkpoints
-*.ipynb
\ No newline at end of file
+*.ipynb
+.devcontainer
+package-lock.json
+package.json
+/node_modules
\ No newline at end of file
diff --git a/docs/make.jl b/docs/make.jl
index d9ac8be98..028d99cac 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -19,7 +19,7 @@ DocMeta.setdocmeta!(QuantumToolbox, :DocTestSetup, doctest_setup; recursive = tr
 
 # some options for `makedocs`
 const DRAFT = get(ENV, "DRAFT", false) == "true"  # `DRAFT   = true`  disables cell evaluation
-const DOCTEST = get(ENV, "DOCTEST", true) == true # `DOCTEST = false` skips doc tests
+const DOCTEST = get(ENV, "DOCTEST", true) == false # `DOCTEST = false` skips doc tests
 
 # generate bibliography
 bib = CitationBibliography(
diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl
index c3d9be754..73165da06 100644
--- a/src/QuantumToolbox.jl
+++ b/src/QuantumToolbox.jl
@@ -100,7 +100,6 @@ include("qobj/operators.jl")
 include("qobj/superoperators.jl")
 include("qobj/synonyms.jl")
 include("qobj/block_diagonal_form.jl")
-
 # time evolution
 include("time_evolution/time_evolution.jl")
 include("time_evolution/callback_helpers/callback_helpers.jl")
@@ -117,6 +116,7 @@ include("time_evolution/mcsolve.jl")
 include("time_evolution/ssesolve.jl")
 include("time_evolution/smesolve.jl")
 include("time_evolution/time_evolution_dynamical.jl")
+include("time_evolution/propagator.jl")
 
 # Others
 include("correlations.jl")
diff --git a/src/time_evolution/propagator.jl b/src/time_evolution/propagator.jl
new file mode 100644
index 000000000..f55f9be12
--- /dev/null
+++ b/src/time_evolution/propagator.jl
@@ -0,0 +1,212 @@
+export Propagator
+export propagator
+
+function propagatorProblem(
+    H::AbstractQuantumObject{Operator},
+    tspan::AbstractVector;
+    params = NullParameters(),
+    progress_bar::Union{Val,Bool} = Val(true),
+    inplace::Union{Val,Bool} = Val(true),
+    kwargs...,
+)
+    H_evo = _sesolve_make_U_QobjEvo(H)
+    U = H_evo.data
+
+    p0 = one(H_evo(0)).data
+
+    kwargs2 = _merge_saveat(tspan, nothing, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...)
+    kwargs3 = _generate_se_me_kwargs(nothing, makeVal(progress_bar), tspan, kwargs2, SaveFuncSESolve)
+
+    return ODEProblem{getVal(inplace),FullSpecialize}(U, p0, tspan, params; kwargs3...)
+end
+
+function propagatorProblem(
+    H::AbstractQuantumObject{SuperOperator},
+    tspan::AbstractVector;
+    params = NullParameters(),
+    progress_bar::Union{Val,Bool} = Val(true),
+    inplace::Union{Val,Bool} = Val(true),
+    kwargs...,
+)
+    L_evo = _mesolve_make_L_QobjEvo(H, [])
+    U = L_evo.data
+
+    p0 = one(L_evo(0)).data
+
+    kwargs2 = _merge_saveat(tspan, nothing, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...)
+    kwargs3 = _generate_se_me_kwargs(nothing, makeVal(progress_bar), tspan, kwargs2, SaveFuncMESolve)
+
+    return ODEProblem{getVal(inplace),FullSpecialize}(U, p0, tspan, params; kwargs3...)
+end
+
+@doc raw"""
+Propagator{T}
+
+Container for time-evolution propagators built from `H` (a Hamiltonian or a Louivillian superoperator).
+This immutable struct groups the generator, time list and solver settings together with any cached
+propagators so that previous results can be reused where applicable.
+
+# Parameters
+- `H::T`
+    The generator of dynamics. `T` is typically an `AbstractQuantumObject{Operator}` for closed,
+    unitary dynamics (Hamiltonian) or an `AbstractQuantumObject{SuperOperator}` for open-system
+    dynamics (Liouvillian or other superoperator).
+- `times::AbstractArray{Float64}`
+    A monotonic array of time points at which propagators are requested or saved. Values are
+    interpreted in the same time units used to build `H`.
+- `props::AbstractArray{A} where A <: AbstractQuantumObject`
+    An array holding computed propagators corresponding to `times`. Elements are quantum objects
+    (operators or superoperators) matching the expected output of propagating with `H`. Both `times` and `props` are initialized with the identity operator/superoperator at t=0.
+- `tol::Float64`
+    This is the time tolerance used to determine if a previous result should be reused. If \$t - t_{\text{cached}} < \text{tol}\$, where $t_\text{cached}\$ is the closed cached
+    time to \$t\$, then the propagator at \$t_\text{cached}\$ will be reused.
+- `solver_kwargs::Base.Pairs`
+    Additional keyword arguments forwarded to the ODE integrator (for example `abstol`, `reltol`,
+    `saveat`, or solver-specific options). Stored as a `Base.Pairs` collection for convenient
+    dispatch into solver APIs.
+- `dims::Union{AbstractArray, Tuple}`
+    Dimension metadata of `H`.
+- `solver::OrdinaryDiffEqAlgorithm`
+    The chosen ODE algorithm (a concrete algorithm type from OrdinaryDiffEq, the default is `Tsit5`) that will be used to integrate the generator to obtain
+    time-ordered evolution.
+- `type`
+    Whether the propagator is an operator or superoperator. 
+
+# Initialization
+A propagator is initialized by calling the `propagator` function.
+
+# Calling
+If `U` is a propagator object, to get the propagator between t0 and t, call `U(t, t0=t0; just_interval::Bool=false)`.
+    - t0 is optional, if left out then t0=0.0. 
+    - `just_interval` determines how the propagator from t0 to t is calculated. 
+        - if `just_interval = false`, then the propagator is calculated by getting U1(t0,0.0) and U2(t, 0.0) and setting
+        ```math
+        U = U2U1^{-1}
+        ```
+        This is useful if U1 is alread cached or one plans to reuse U1 or U2. 
+        - if `just_interval=true` then the propagator is calulated by directly integrating from t0 to t. This approach does't save 
+        anything but is useful if t-t0 >> t0 as it avoids the long integration required to get t0. 
+"""
+struct Propagator{T<:Union{AbstractQuantumObject{Operator},AbstractQuantumObject{SuperOperator}}}
+    H::T
+    times::AbstractArray{Float64}
+    props::AbstractArray{A} where {A<:AbstractQuantumObject}
+    tol::Float64
+    solver_kwargs::Base.Pairs
+    dims::Union{AbstractArray,Tuple}
+    solver::OrdinaryDiffEqAlgorithm
+    type::Any
+end
+
+function Base.show(io::IO, p::Propagator)
+    println("Propagator: ")
+    println("   dims: $(p.dims)")
+    println("   issuper: $(issuper(p.H))")
+    println("   size: $(size(p.H))")
+    return println("   Number of Saved Times: $(length(p.times))")
+end
+
+@doc raw"""
+propagator(H::Union{QobjEvo, QuantumObject}; tol=1e-6, solver=Tsit5(), kwargs...)
+
+Construct and return a Propagator object prepared for time evolution using the
+given Hamiltonian or time-dependent operator.
+
+# Arguments
+- H::Union{QobjEvo, QuantumObject}
+    The Hamiltonian or generator for the dynamics. Can be a time-independent
+    QuantumObject or a time-dependent QobjEvo describing H(t).
+- c_ops::Union{AbstractArray, Tuple} (Optional)
+    - if no c_ops are given, then this is treated as a lossless propagator and the propagator itself will be of type `AbstractQuantumObject{Operator}`. If c_ops are given, 
+    then the propagator will be of type `AbstractQuantumObject{SuperOperator}` and the save `H` will be the Louivillian.
+- tol::Real=1e-6
+    Absolute/relative tolerance used by the ODE integrator. This value is
+    forwarded to the integrator/Propagator to control numerical accuracy.
+- solver
+    The ODE solver algorithm to use (default: `Tsit5()`); any solver object
+    compatible with the underlying ODE interface is accepted.
+- kwargs...
+    Additional keyword arguments are forwarded to the Propagator constructor
+    (e.g., options for step control, callback functions, or integration
+    settings supported by the backend).
+
+# Notes
+- The function does not perform propagation itself; it only constructs and
+  configures the Propagator object which can then be used to evolve states or
+  operators.
+- For time-dependent `H` provide a `QobjEvo` with an appropriate time
+  dependence. For time-independent dynamics provide a `QuantumObject`.
+
+# Returns
+- `Propagator`
+    An initialized Propagator instance set up with
+    - initial time [0.0],
+    - initial operator equal to the identity on the Hilbert space of `H`
+      (created with `one(H)`),
+    - the requested tolerance and solver,
+    - dimensions inferred from `H`, and
+    - any extra options passed via `kwargs`.
+"""
+function propagator(H::Union{QobjEvo,QuantumObject}; tol = 1e-6, solver = Tsit5(), kwargs...)
+    H_evo = QobjEvo(H)
+    return Propagator(H, [0.0], [one(H_evo(0))], tol, kwargs, H.dims, solver, Operator())
+end
+
+function propagator(
+    H::Union{QobjEvo,QuantumObject},
+    c_ops::Union{AbstractArray,Tuple};
+    tol = 1e-6,
+    solver = Tsit5(),
+    kwargs...,
+)
+    L = liouvillian(H, c_ops)
+    L_evo = QobjEvo(L)
+    return Propagator(L, [0.0], [one(L_evo(0))], tol, kwargs, L.dims, solver, SuperOperator())
+end
+
+function _lookup_or_compute(p::Propagator, t::Real)
+    """
+    Get U(t) from cache or compute it.
+    """
+    idx = searchsorted(p.times, t)
+
+    # Check if we found an exact match or close enough within tolerance
+    if (idx.start <= length(p.times)) && (abs(t - p.times[idx.start]) <= p.tol)
+        U = p.props[idx.start]
+    elseif (idx.start > 1) && (abs(t - p.times[idx.start-1]) <= p.tol)
+        U = p.props[idx.start-1]
+    else
+        t0 = (idx.start > 1) ? p.times[idx.start-1] : 0.0
+        U = _compute_propagator(p, t)
+        insert!(p.props, idx.start, U)
+        insert!(p.times, idx.start, t)
+    end
+
+    return U
+end
+
+function _compute_propagator(p::Propagator, t::Real; t0 = 0.0)
+    prob = propagatorProblem(p.H, [t0, t]; p.solver_kwargs...)
+    res = solve(prob, p.solver)
+    U = QuantumObject(res.u[end], dims = p.dims, type = p.type)
+    return U
+end
+
+function (p::Propagator)(t::Real, t0 = 0.0; just_interval = false)
+    if just_interval
+        return _compute_propagator(p, t; t0 = t0)
+    end
+    # We start by computing U_t0 if needed. This means that the computation for U will, at worst, start at t0 instead of 0.0.
+    if t0 != 0.0
+        U_t0 = (_lookup_or_compute(p, t0))
+    end
+    U = _lookup_or_compute(p, t)
+
+    # We only do the multiplication if needed
+    if t0 != 0.0
+        U = U*inv(U_t0)
+    end
+    return U
+end
+
+Base.size(p::Propagator) = size(p.H)
diff --git a/test/core-test/propagator.jl b/test/core-test/propagator.jl
new file mode 100644
index 000000000..1f07b7ba8
--- /dev/null
+++ b/test/core-test/propagator.jl
@@ -0,0 +1,48 @@
+@testitem "Propagators" begin
+    sx = sigmax()
+    sy = sigmay()
+    sz = sigmaz()
+    I = one(sz)
+
+    H = sz + QobjEvo((0.1*sx, (p, t) -> sin(2*t)))
+
+    p_se = propagator(H)
+
+    c_ops = [0.1*sz]
+    p_me = propagator(H, c_ops)
+
+    sesolve_res = []
+    mesolve_res = []
+
+    t0 = 0.0
+    t = 10.0
+    for i in 0:1
+        psi = fock(2, i)
+        push!(sesolve_res, sesolve(H, psi, [t0, t]))
+    end
+
+    sup_psi_base = mat2vec(fock(2, 0)*fock(2, 0)')
+    for i in 1:4
+        sup_psi = 0 * sup_psi_base
+        sup_psi[i] = 1
+        push!(mesolve_res, mesolve(H, vec2mat(sup_psi), [t0, t], c_ops))
+    end
+
+    U_se = p_se(t)
+    U_me = p_me(t)
+
+    U_sesolve = QuantumObject(hcat([sesolve_res[i].states[end].data for i in 1:2]...))
+    U_mesolve = QuantumObject(hcat([mat2vec(mesolve_res[i].states[end]).data for i in 1:4]...), type = SuperOperator())
+
+    @test norm(U_se-U_sesolve) < 1e-8
+    @test norm(U_me-U_mesolve) < 1e-8
+
+    U_se = p_se(t, 3.0)
+    U_se_JI = p_se(t, 3.0; just_interval = true)
+    U_me = p_me(t, 3.0)
+    U_me_JI = p_me(t, 3.0; just_interval = true)
+
+    @test norm(U_se-U_se_JI)<1e-6
+    @test norm(U_me-U_me_JI) < 1e-6
+
+end