diff --git a/Project.toml b/Project.toml
index 0f73682..4bc32fc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "PiccoloQuantumObjects"
 uuid = "5a402ddf-f93c-42eb-975e-5582dcda653d"
 authors = ["Aaron Trowbridge <aaron.j.trowbridge@gmail.com> and contributors"]
-version = "0.6.1"
+version = "0.7.0"
 
 [deps]
 ExponentialAction = "e24c0720-ea99-47e8-929e-571b494574d3"
diff --git a/docs/Project.toml b/docs/Project.toml
index 4c09eed..0e91532 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
@@ -8,4 +9,10 @@ PiccoloQuantumObjects = "5a402ddf-f93c-42eb-975e-5582dcda653d"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 
 [compat]
-Literate = "2.7.0"
+CairoMakie = "0.13"
+Documenter = "1.15"
+Literate = "2.20"
+LiveServer = "1.5"
+NamedTrajectories = "0.5"
+PiccoloDocsTemplate = "0.2"
+Revise = "3.11"
diff --git a/docs/literate/quantum_systems.jl b/docs/literate/quantum_systems.jl
index 5a1407b..e897081 100644
--- a/docs/literate/quantum_systems.jl
+++ b/docs/literate/quantum_systems.jl
@@ -22,24 +22,30 @@ The [`QuantumSystem`](@ref) type is used to represent a quantum system with a dr
 Hamiltonian and a set of drive Hamiltonians,
 
 ```math
-H = H_{\text{drift}} + \sum_i a_i H_{\text{drives}}^{(i)}
+H(u, t) = H_{\text{drift}} + \sum_i u_i H_{\text{drives}}^{(i)}
 ```
 
+where ``u`` is the control vector and ``t`` is time.
+
 ```@docs; canonical = false
 QuantumSystem
 ```
 
 `QuantumSystem`'s are containers for quantum dynamics. Internally, they compute the
-necessary isomorphisms to perform the dynamics in a real vector space.
+necessary isomorphisms to perform the dynamics in a real vector space. All systems
+require explicit specification of `T_max` (maximum time) and `drive_bounds` (control bounds).
 
 =#
 
 H_drift = PAULIS[:Z]
 H_drives = [PAULIS[:X], PAULIS[:Y]]
-system = QuantumSystem(H_drift, H_drives)
+T_max = 10.0
+drive_bounds = [(-1.0, 1.0), (-1.0, 1.0)]
+system = QuantumSystem(H_drift, H_drives, T_max, drive_bounds)
 
-a_drives = [1, 0]
-system.H(a_drives)
+u_controls = [1.0, 0.0]
+t = 0.0
+system.H(u_controls, t)
 
 #=
 To extract the drift and drive Hamiltonians from a `QuantumSystem`, use the 
@@ -58,27 +64,13 @@ drives[2] |> sparse
 
 #=
 !!! note
-    We can also construct a `QuantumSystem` directly from a Hamiltonian function. Internally,
-    `ForwardDiff.jl` is used to compute the drives.
+    We can also construct a `QuantumSystem` directly from a Hamiltonian function.
+    The function must accept `(u, t)` arguments where `u` is the control vector and `t` is time.
 =#
 
-H(a) = PAULIS[:Z] + a[1] * PAULIS[:X] + a[2] * PAULIS[:Y]
-system = QuantumSystem(H, 2)
-get_drives(system)[1] |> sparse
-
-# _Create a noise model with a confusion matrix._
-function H(a; C::Matrix{Float64}=[1.0 0.0; 0.0 1.0])
-    b = C * a
-    return b[1] * PAULIS.X + b[2] * PAULIS.Y
-end
-
-C_matrix = [0.99 0.01; -0.01 1.01]
-system = QuantumSystem(a -> H(a, C=C_matrix), 2; params=Dict(:C => C_matrix))
-confused_drives = get_drives(system)
-confused_drives[1] |> sparse
-
-# 
-confused_drives[2] |> sparse
+H(u, t) = PAULIS[:Z] + u[1] * PAULIS[:X] + u[2] * PAULIS[:Y]
+system = QuantumSystem(H, 10.0, [(-1.0, 1.0), (-1.0, 1.0)])
+system.H([1.0, 0.0], 0.0) |> sparse
 
 #=
 ## Open quantum systems
@@ -95,7 +87,9 @@ OpenQuantumSystem
 H_drives = [PAULIS[:X]]
 a = annihilate(2)
 dissipation_operators = [a'a, a]
-system = OpenQuantumSystem(H_drives, dissipation_operators=dissipation_operators)
+T_max = 10.0
+drive_bounds = [(-1.0, 1.0)]
+system = OpenQuantumSystem(H_drives, T_max, drive_bounds, dissipation_operators=dissipation_operators)
 system.dissipation_operators[1] |> sparse
 
 # 
@@ -110,41 +104,6 @@ system.dissipation_operators[2] |> sparse
 
 get_drift(system) |> sparse
 
-#=
-## Time Dependent Quantum Systems
-A [`TimeDependentQuantumSystem`](@ref) is a `QuantumSystem` with time-dependent Hamiltonians.
-```@docs; canonical = false
-TimeDependentQuantumSystem
-```
-
-A function `H(a, t)` or carrier and phase kwargs are used to specify time-dependent drives,
-```math
-    H(a, t) = H_{\text{drift}} + \sum_i a_i \cos(\omega_i t + \phi_i) H_{\text{drives}}^{(i)}
-```
-=#
-# _Create a time-dependent Hamiltonian with a time-dependent drive._
-H(a, t) = PAULIS.Z + a[1] * cos(t) * PAULIS.X
-system = TimeDependentQuantumSystem(H, 1)
-
-# _The drift Hamiltonian is the Z operator, but its now a function of time!_
-get_drift(system)(0.0) |> sparse
-
-# _The drive Hamiltonian is the X operator, but its now a function of time!_
-get_drives(system)[1](0.0) |> sparse
-
-# _Change the time to π._
-get_drives(system)[1](π) |> sparse
-
-# _Similar matrix constructors exist, but with carrier and phase kwargs._
-system = TimeDependentQuantumSystem(PAULIS.Z, [PAULIS.X], carriers=[1.0], phases=[0.0])
-
-# _This is the same as before, t=0.0:_
-get_drives(system)[1](0.0) |> sparse
-
-# _and at π:_
-get_drives(system)[1](π) |> sparse
-
-
 #=
 ## Composite quantum systems
 
@@ -155,10 +114,10 @@ lifted to the full Hilbert space.
 
 =#
 
-system_1 = QuantumSystem([PAULIS[:X]])
-system_2 = QuantumSystem([PAULIS[:Y]])
+system_1 = QuantumSystem([PAULIS[:X]], 1.0, [(-1.0, 1.0)])
+system_2 = QuantumSystem([PAULIS[:Y]], 1.0, [(-1.0, 1.0)])
 H_drift = PAULIS[:Z] ⊗ PAULIS[:Z]
-system = CompositeQuantumSystem(H_drift, [system_1, system_2]);
+system = CompositeQuantumSystem(H_drift, Matrix{ComplexF64}[], [system_1, system_2], 1.0, Float64[]);
 
 # _The drift Hamiltonian is the ZZ coupling._
 get_drift(system) |> sparse
@@ -214,11 +173,11 @@ is_reachable
 =#
 
 # _Y can be reached by commuting Z and X._
-system = QuantumSystem(PAULIS[:Z], [PAULIS[:X]])
+system = QuantumSystem(PAULIS[:Z], [PAULIS[:X]], 1.0, [(-1.0, 1.0)])
 is_reachable(PAULIS[:Y], system)
 
 # _Y cannot be reached by X alone._
-system = QuantumSystem([PAULIS[:X]])
+system = QuantumSystem([PAULIS[:X]], 1.0, [(-1.0, 1.0)])
 is_reachable(PAULIS[:Y], system)
 
 #=
@@ -231,7 +190,7 @@ direct_sum
 =#
 
 # _Create a pair of non-interacting qubits._
-system_1 = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])
-system_2 = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])
+system_1 = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
+system_2 = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
 system = direct_sum(system_1, system_2)
 get_drift(system) |> sparse
diff --git a/docs/literate/quickstart.jl b/docs/literate/quickstart.jl
new file mode 100644
index 0000000..640eff2
--- /dev/null
+++ b/docs/literate/quickstart.jl
@@ -0,0 +1,143 @@
+# # Quickstart
+
+# This quickstart guide will help you get up and running with PiccoloQuantumObjects.jl.
+
+# ## Installation
+
+# ```julia
+# using Pkg
+# Pkg.add("PiccoloQuantumObjects")
+# ```
+
+# ## Basic Usage
+
+using PiccoloQuantumObjects
+using LinearAlgebra
+using SparseArrays
+using NamedTrajectories
+
+# ### Creating a Quantum System
+
+# Define Hamiltonian components
+H_drift = PAULIS[:Z]
+H_drives = [PAULIS[:X], PAULIS[:Y]]
+
+# Specify time and control bounds
+T_max = 10.0  # Maximum evolution time
+drive_bounds = [(-1.0, 1.0), (-1.0, 1.0)]  # Control bounds for each drive
+
+# Create the quantum system
+system = QuantumSystem(H_drift, H_drives, T_max, drive_bounds)
+
+# ### Working with Quantum States
+
+# Create quantum states
+ψ_ground = ket_from_string("g", [2])  # Ground state
+
+#
+ψ_excited = ket_from_string("e", [2])  # Excited state
+
+# Calculate fidelity
+fidelity(ψ_ground, ψ_excited)
+
+# ### Simulating Evolution (Rollouts)
+
+# Initial state
+ψ_init = ComplexF64[1.0, 0.0]
+
+# Define control sequence
+T = 10  # Number of time steps
+controls = rand(2, T)  # Random controls for two drives
+Δt = fill(0.1, T)  # Time step duration
+
+# Perform rollout
+ψ̃_rollout = rollout(ψ_init, controls, Δt, system)
+
+# Check final state fidelity
+ψ_goal = ComplexF64[0.0, 1.0]
+rollout_fidelity(ψ_init, ψ_goal, controls, Δt, system)
+
+# ### Open Quantum Systems
+
+# Add dissipation operators
+a = annihilate(2)
+dissipation_operators = [a'a, a]
+
+# Create open quantum system
+open_system = OpenQuantumSystem(
+    H_drives, 
+    T_max, 
+    drive_bounds, 
+    dissipation_operators=dissipation_operators
+)
+
+# ### Composite Systems
+
+# Create subsystems
+sys1 = QuantumSystem([PAULIS[:X]], 10.0, [(-1.0, 1.0)])
+sys2 = QuantumSystem([PAULIS[:Y]], 10.0, [(-1.0, 1.0)])
+
+# Define coupling
+H_coupling = 0.1 * kron(PAULIS[:Z], PAULIS[:Z])
+
+# Create composite system
+composite_sys = CompositeQuantumSystem(
+    H_coupling, 
+    Matrix{ComplexF64}[], 
+    [sys1, sys2], 
+    10.0, 
+    Float64[]
+)
+
+# ## Visualization
+
+# PiccoloQuantumObjects.jl integrates with NamedTrajectories.jl for plotting trajectories.
+
+using CairoMakie
+
+# ### Plotting Controls and States
+
+# Create a trajectory with controls and states
+T_plot = 50
+controls_plot = 0.5 * sin.(2π * (1:T_plot) / T_plot)
+Δt_plot = fill(T_max / T_plot, T_plot)
+
+# Perform a rollout to get the state evolution
+ψ̃_traj = rollout(ψ_init, hcat(controls_plot, -controls_plot)', Δt_plot, system)
+
+# Create a NamedTrajectory for plotting
+traj = NamedTrajectory(
+    (
+        ψ̃ = ψ̃_traj,
+        a = hcat(controls_plot, -controls_plot)',
+        Δt = Δt_plot
+    );
+    timestep=:Δt,
+    controls=:a,
+    initial=(ψ̃ = ket_to_iso(ψ_init),),
+    goal=(ψ̃ = ket_to_iso(ψ_goal),)
+)
+
+# Plot the trajectory 
+plot(traj)
+
+# ### Plotting State Populations
+
+# Use transformations to plot populations directly from isomorphic states
+plot(
+    traj, 
+    [:a],
+    transformations=[
+        :ψ̃ => (ψ̃ -> abs2.(iso_to_ket(ψ̃)))
+    ],
+    transformation_labels=["Populations"]
+)
+
+#=
+## Next Steps
+
+- Explore the [Quantum Systems](@ref) manual for detailed system construction
+- Learn about [Quantum Objects](@ref) for working with states and operators
+- See [Rollouts](@ref) for simulation and fidelity calculations
+- Understand [Isomorphisms](@ref) for the underlying mathematical transformations
+=#
diff --git a/docs/literate/rollouts.jl b/docs/literate/rollouts.jl
index 5364ae4..b488db9 100644
--- a/docs/literate/rollouts.jl
+++ b/docs/literate/rollouts.jl
@@ -81,7 +81,7 @@ T = 10
 ψ_init = ComplexF64[1.0, 0.0]
 controls = rand(2, T)
 Δt = fill(0.1, T)
-system = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])
+system = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
 ψ̃_rollout = rollout(ψ_init, controls, Δt, system) 
 ψ̃_rollout |> size
 
diff --git a/docs/make.jl b/docs/make.jl
index b166af9..5f2b71a 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -3,11 +3,12 @@ using PiccoloDocsTemplate
 
 pages = [
     "Home" => "index.md",
+    "Quickstart" => "generated/quickstart.md",
     "Manual" => [
-        "Isomorphisms" => "generated/isomorphisms.md",
-        "Quantum Objects" => "generated/quantum_objects.md",
         "Quantum Systems" => "generated/quantum_systems.md",
+        "Quantum Objects" => "generated/quantum_objects.md",
         "Rollouts" => "generated/rollouts.md",
+        "Isomorphisms" => "generated/isomorphisms.md",
     ],
     "Library" => "lib.md",
 ]
diff --git a/ext/QuantumToolboxExt.jl b/ext/QuantumToolboxExt.jl
index 8338d4b..6df84db 100644
--- a/ext/QuantumToolboxExt.jl
+++ b/ext/QuantumToolboxExt.jl
@@ -71,10 +71,10 @@ get_c_ops(sys::OpenQuantumSystem) = Qobj.(sys.dissipation_operators)
     # X gate
     traj = named_trajectory_type_2()
 
-    sys = QuantumSystem([PAULIS.X, PAULIS.Y])
+    sys = QuantumSystem([PAULIS.X, PAULIS.Y], 3.92, [(-1.0, 1.0), (-1.0, 1.0)])
 
     open_sys = OpenQuantumSystem(
-        [PAULIS.X, PAULIS.Y], dissipation_operators=[annihilate(2)]
+        [PAULIS.X, PAULIS.Y], 1.0, [1.0, 1.0], dissipation_operators=[annihilate(2)]
     )
 
     ψ0 = Qobj([1; 0])
diff --git a/src/direct_sums.jl b/src/direct_sums.jl
index dc5ca21..6a6f8bf 100644
--- a/src/direct_sums.jl
+++ b/src/direct_sums.jl
@@ -3,9 +3,11 @@ module DirectSums
 export direct_sum
 
 using SparseArrays
+using SparseArrays: blockdiag
 using TestItems
 
 using ..QuantumSystems
+using ..Isomorphisms: operator_to_iso_vec, iso_vec_to_operator
 
 
 """
@@ -44,36 +46,30 @@ end
     direct_sum(sys1::QuantumSystem, sys2::QuantumSystem)
 
 Returns the direct sum of two `QuantumSystem` objects.
+
+Constructs a new system where the Hilbert space is the direct sum of the two input systems:
+H = H₁ ⊕ H₂ = [H₁  0 ]
+               [0   H₂]
+
+Both systems must have the same number of drives. The resulting system uses sys1's T_max and drive_bounds.
+
+# Example
+```julia
+sys1 = QuantumSystem([PAULIS[:X]], 10.0, [(-1.0, 1.0)])
+sys2 = QuantumSystem([PAULIS[:Y]], 10.0, [(-1.0, 1.0)])
+sys_combined = direct_sum(sys1, sys2)
+```
 """
 function direct_sum(sys1::QuantumSystem, sys2::QuantumSystem)
     @assert sys1.n_drives == sys2.n_drives "System 1 drives ($(sys1.n_drives)) must equal System 2 drives ($(sys2.n_drives))"
     n_drives = sys1.n_drives
-    H = a -> direct_sum(sys1.H(a), sys2.H(a))
-    direct_sum_params = Dict{Symbol, Dict{Symbol, Any}}()
-    if haskey(sys1.params, :system_1)
-        n_systems = length(keys(sys1.params))
-        direct_sum_params = sys1.params
-        if haskey(sys2.params, :system_1)
-            for i = 1:length(keys(sys2.params))
-                direct_sum_params[Symbol("system_$(n_systems + i)")] =
-                    sys2.params[Symbol("system_$(i)")]
-            end
-        else
-            direct_sum_params[Symbol("system_$(n_systems + 1)")] = sys2.params
-        end
-    else
-        direct_sum_params[:system_1] = sys1.params
-        if haskey(sys2.params, :system_1)
-            n_systems = length(keys(sys2.params))
-            for i = 1:length(keys(sys2.params))
-                direct_sum_params[Symbol("system_$(1 + i)")] =
-                    sys2.params[Symbol("system_$(i)")]
-            end
-        else
-            direct_sum_params[:system_2] = sys2.params
-        end
-    end
-    return QuantumSystem(H, n_drives; params=direct_sum_params)
+    H = (u, t) -> direct_sum(sys1.H(u, t), sys2.H(u, t))
+    
+    # Get combined drive bounds from both systems - but they should be the same length since n_drives must match
+    drive_bounds = sys1.drive_bounds  # They should be the same as sys2.drive_bounds
+    
+    # Create new QuantumSystem with the Hamiltonian function
+    return QuantumSystem(H, sys1.T_max, drive_bounds)
 end
 
 direct_sum(systems::AbstractVector{<:QuantumSystem}) = reduce(direct_sum, systems)
@@ -93,21 +89,11 @@ direct_sum(systems::AbstractVector{<:QuantumSystem}) = reduce(direct_sum, system
 end
 
 @testitem "Test quantum system direct sum" begin
-    sys1 = QuantumSystem([1 2; 3 4])
-    sys2 = QuantumSystem([5 6; 7 8])
-    sys = direct_sum(sys1, sys2)
-    @test sys.H(0) == [1 2 0 0; 3 4 0 0; 0 0 5 6; 0 0 7 8]
-    @test sys.n_drives == 0
-    # created new keys for the direct sum
-    @test sys.params[:system_1] == Dict()
-    @test sys.params[:system_2] == Dict()
-
-    sys1 = QuantumSystem([1 2; 3 4]; params=Dict(:a => 1))
-    sys2 = QuantumSystem([5 6; 7 8]; params=Dict(:b => 2))
+    sys1 = QuantumSystem(ComplexF64[1 2; 3 4], 1.0)
+    sys2 = QuantumSystem(ComplexF64[5 6; 7 8], 1.0)
     sys = direct_sum(sys1, sys2)
-    @test sys.H(0) == [1 2 0 0; 3 4 0 0; 0 0 5 6; 0 0 7 8]
+    @test sys.H(Float64[], 0.0) == ComplexF64[1 2 0 0; 3 4 0 0; 0 0 5 6; 0 0 7 8]
     @test sys.n_drives == 0
-    @test sys.params == Dict{Symbol, Dict{Symbol, Any}}(:system_1 => Dict(:a => 1), :system_2 => Dict(:b => 2))
 end
 
 end
diff --git a/src/embedded_operators.jl b/src/embedded_operators.jl
index bbf15c1..3f9d78f 100644
--- a/src/embedded_operators.jl
+++ b/src/embedded_operators.jl
@@ -491,10 +491,12 @@ end
 
 @testitem "Embedded operator from system" begin
     CZ = GATES[:CZ]
-    a = [0 1 0; 0 0 1; 0 0 0]
+    a = annihilate(3) 
     σ_x = a + a'
     σ_y = -1im*(a - a')
-    system = QuantumSystem([kron(σ_x, σ_x), kron(σ_y, σ_y)])
+    T_max = 1.0
+    u_bounds = ones(2)
+    system = QuantumSystem([kron(σ_x, σ_x), kron(σ_y, σ_y)], T_max, u_bounds)
 
     op_explicit_qubit = EmbeddedOperator(
         CZ,
@@ -629,7 +631,7 @@ end
     @test embedded_op.operator == kron(embed(I(2), 1:2, 3), embed(subspace_op, subspace_op_indices, 3^2))
 
     # Composite system
-    system = CompositeQuantumSystem([QuantumSystem([P]) for P ∈ PAULIS])
+    system = CompositeQuantumSystem([QuantumSystem([P], 1.0, [(-1.0, 1.0)]) for P ∈ PAULIS], 1.0, Float64[])
     embedded_op = EmbeddedOperator(subspace_op, [2, 3], fill(1:2, length(PAULIS)), system)
     # 4 PAULIS
     @test embedded_op.operator == kron(I(2), subspace_op, I(2))
diff --git a/src/quantum_system_utils.jl b/src/quantum_system_utils.jl
index 1c60e08..47ae396 100644
--- a/src/quantum_system_utils.jl
+++ b/src/quantum_system_utils.jl
@@ -336,12 +336,12 @@ end
     @test is_reachable(target, gen, compute_basis=true, verbose=false)
 
     # System
-    sys = QuantumSystem([PAULIS[:X], PAULIS[:Y], PAULIS[:Z]])
+    sys = QuantumSystem([PAULIS[:X], PAULIS[:Y], PAULIS[:Z]], 1.0, [(-1.0, 1.0), (-1.0, 1.0), (-1.0, 1.0)])
     target = PAULIS[:Z]
     @test is_reachable(target, sys, verbose=false)
 
     # System with drift
-    sys = QuantumSystem(PAULIS[:Z], [PAULIS[:X]])
+    sys = QuantumSystem(PAULIS[:Z], [PAULIS[:X]], 1.0, [(-1.0, 1.0)])
     target = PAULIS[:Z]
     @test is_reachable(target, sys, verbose=false)
 end
@@ -382,8 +382,8 @@ end
     @test is_reachable(XI, incomplete_gen, verbose=false)
 
     # QuantumSystems
-    complete_gen_sys = QuantumSystem(complete_gen)
-    incomplete_gen_sys = QuantumSystem(incomplete_gen)
+    complete_gen_sys = QuantumSystem(complete_gen, 1.0, fill((-1.0, 1.0), length(complete_gen)))
+    incomplete_gen_sys = QuantumSystem(incomplete_gen, 1.0, fill((-1.0, 1.0), length(incomplete_gen)))
     # Pass
     @test is_reachable(R2, complete_gen_sys, verbose=false)
     @test is_reachable(CZ, complete_gen_sys, verbose=false)
@@ -404,7 +404,7 @@ end
     @test is_reachable(target, gen, verbose=false)
 
     # System
-    sys = QuantumSystem([GATES[:X], GATES[:Y]])
+    sys = QuantumSystem([GATES[:X], GATES[:Y]], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
     @test is_reachable(target, sys, verbose=false)
 end
 
diff --git a/src/quantum_systems/_quantum_systems.jl b/src/quantum_systems/_quantum_systems.jl
index 30027c6..707346f 100644
--- a/src/quantum_systems/_quantum_systems.jl
+++ b/src/quantum_systems/_quantum_systems.jl
@@ -34,7 +34,7 @@ abstract type AbstractQuantumSystem end
 
 Returns the drift Hamiltonian of the system.
 """
-get_drift(sys::AbstractQuantumSystem) = sys.H(zeros(sys.n_drives))
+get_drift(sys::AbstractQuantumSystem) = sys.H(zeros(sys.n_drives), 0.0)
 
 """
     get_drives(sys::AbstractQuantumSystem)
@@ -44,7 +44,7 @@ Returns the drive Hamiltonians of the system.
 function get_drives(sys::AbstractQuantumSystem)
     H_drift = get_drift(sys)
     # Basis vectors for controls will extract drive operators
-    return [sys.H(I[1:sys.n_drives, i]) - H_drift for i ∈ 1:sys.n_drives]
+    return [sys.H(I[1:sys.n_drives, i], 0.0) - H_drift for i ∈ 1:sys.n_drives]
 end
 
 function Base.show(io::IO, sys::AbstractQuantumSystem)
@@ -63,6 +63,6 @@ function get_c_ops end
 
 include("quantum_systems.jl")
 include("composite_quantum_systems.jl")
-include("time_dependent_quantum_systems.jl")
+# include("time_dependent_quantum_systems.jl")  # Commented out - not yet refactored
 
 end
diff --git a/src/quantum_systems/composite_quantum_systems.jl b/src/quantum_systems/composite_quantum_systems.jl
index ceaa767..439f00d 100644
--- a/src/quantum_systems/composite_quantum_systems.jl
+++ b/src/quantum_systems/composite_quantum_systems.jl
@@ -95,103 +95,232 @@ end
 """
     CompositeQuantumSystem <: AbstractQuantumSystem
 
-A composite quantum system consisting of `subsystems`. Couplings between subsystems can
-be additionally defined. Subsystem drives are always appended to any new coupling drives.
+A composite quantum system consisting of multiple `subsystems` with optional coupling terms.
+
+Composite systems represent multiple quantum subsystems (e.g., multiple qubits or oscillators)
+that may be coupled together. Each subsystem's Hamiltonians are automatically lifted to the 
+full tensor product space, and subsystem drives are appended to any coupling drives.
+
+# Fields
+- `H::Function`: The total Hamiltonian function: (u, t) -> H(u, t)
+- `G::Function`: The isomorphic generator function: (u, t) -> G(u, t)
+- `H_drift::SparseMatrixCSC{ComplexF64, Int}`: The total drift Hamiltonian including subsystem drifts and couplings
+- `H_drives::Vector{SparseMatrixCSC{ComplexF64, Int}}`: All drive Hamiltonians (coupling drives + subsystem drives)
+- `T_max::Float64`: Maximum evolution time
+- `drive_bounds::Vector{Tuple{Float64, Float64}}`: Drive amplitude bounds for each control
+- `n_drives::Int`: Total number of control drives
+- `levels::Int`: Total dimension of the composite system (product of subsystem dimensions)
+- `subsystem_levels::Vector{Int}`: Dimensions of each subsystem
+- `subsystems::Vector{QuantumSystem}`: The individual quantum subsystems
+
+See also [`QuantumSystem`](@ref), [`lift_operator`](@ref).
+
+# Example
+```julia
+# Two qubits with ZZ coupling
+sys1 = QuantumSystem([PAULIS[:X]], 10.0, [(-1.0, 1.0)])
+sys2 = QuantumSystem([PAULIS[:Y]], 10.0, [(-1.0, 1.0)])
+H_coupling = 0.1 * kron(PAULIS[:Z], PAULIS[:Z])
+csys = CompositeQuantumSystem(H_coupling, [sys1, sys2], 10.0, Float64[])
+```
 """
 struct CompositeQuantumSystem{F1<:Function, F2<:Function} <: AbstractQuantumSystem
     H::F1
     G::F2
+    H_drift::SparseMatrixCSC{ComplexF64, Int}
+    H_drives::Vector{SparseMatrixCSC{ComplexF64, Int}}
+    T_max::Float64
+    drive_bounds::Vector{Tuple{Float64, Float64}}
     n_drives::Int
     levels::Int
-    params::Dict{Symbol, Any}
     subsystem_levels::Vector{Int}
     subsystems::Vector{QuantumSystem}
+end
 
-    function CompositeQuantumSystem(
-        H_drift::AbstractMatrix{<:Number},
-        H_drives::AbstractVector{<:AbstractMatrix{<:Number}},
-        subsystems::AbstractVector{<:QuantumSystem};
-        params::Dict{Symbol, Any}=Dict{Symbol, Any}()
+"""
+    CompositeQuantumSystem(
+        H_drift::AbstractMatrix,
+        H_drives::AbstractVector{<:AbstractMatrix},
+        subsystems::AbstractVector{<:QuantumSystem},
+        T_max::Float64,
+        drive_bounds::Vector
     )
-        subsystem_levels = [sys.levels for sys ∈ subsystems]
-        levels = prod(subsystem_levels)
 
-        H_drift = sparse(H_drift)
-        for (i, sys) ∈ enumerate(subsystems)
-            H_drift += lift_operator(get_drift(sys), i, subsystem_levels)
-        end
+Construct a CompositeQuantumSystem with coupling drift and drive terms.
+
+# Arguments
+- `H_drift::AbstractMatrix`: Coupling drift Hamiltonian (in full tensor product space)
+- `H_drives::AbstractVector{<:AbstractMatrix}`: Coupling drive Hamiltonians
+- `subsystems::AbstractVector{<:QuantumSystem}`: Vector of subsystems to compose
+- `T_max::Float64`: Maximum evolution time
+- `drive_bounds::Vector`: Drive bounds for the coupling drives (subsystem bounds are inherited)
+
+The total drift includes both the coupling drift and all subsystem drifts (automatically lifted).
+The total drives include coupling drives followed by all subsystem drives (automatically lifted).
+
+# Example
+```julia
+sys1 = QuantumSystem(PAULIS[:Z], [PAULIS[:X]], 10.0, [(-1.0, 1.0)])
+sys2 = QuantumSystem([PAULIS[:Y]], 10.0, [(-1.0, 1.0)])
+g12 = 0.1 * kron(PAULIS[:X], PAULIS[:X])  # coupling drift
+csys = CompositeQuantumSystem(g12, Matrix{ComplexF64}[], [sys1, sys2], 10.0, Float64[])
+```
+"""
+function CompositeQuantumSystem(
+    H_drift::AbstractMatrix{<:Number},
+    H_drives::AbstractVector{<:AbstractMatrix{<:Number}},
+    subsystems::AbstractVector{<:QuantumSystem},
+    T_max::Float64,
+    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}
+)
+    # Normalize drive bounds to tuples
+    drive_bounds = [
+        b isa Tuple ? b : (-b, b) for b in drive_bounds
+    ]
+    
+    subsystem_levels = [sys.levels for sys ∈ subsystems]
+    levels = prod(subsystem_levels)
+
+    H_drift = sparse(H_drift)
+    for (i, sys) ∈ enumerate(subsystems)
+        H_drift += lift_operator(get_drift(sys), i, subsystem_levels)
+    end
 
-        H_drives = sparse.(H_drives)
-        for (i, sys) ∈ enumerate(subsystems)
-            for H_drive ∈ get_drives(sys)
-                push!(H_drives, lift_operator(H_drive, i, subsystem_levels))
-            end
+    H_drives = sparse.(H_drives)
+    for (i, sys) ∈ enumerate(subsystems)
+        for H_drive ∈ get_drives(sys)
+            push!(H_drives, lift_operator(H_drive, i, subsystem_levels))
         end
+    end
 
-        n_drives = length(H_drives)
-        H_drives = sparse.(H_drives)
-        G_drives = sparse.(Isomorphisms.G.(H_drives))
+    n_drives = length(H_drives)
+    H_drives = sparse.(H_drives)
+    G_drift = sparse(Isomorphisms.G(H_drift))
+    G_drives = sparse.(Isomorphisms.G.(H_drives))
+
+    if n_drives == 0
+        H = (u, t) -> H_drift
+        G = (u, t) -> G_drift
+    else
+        H = (u, t) -> H_drift + sum(u .* H_drives)
+        G = (u, t) -> G_drift + sum(u .* G_drives)
+    end
 
-        if n_drives == 0
-            H = a -> H_drift
-            G = a -> Isomorphisms.G(H_drift)
-        else
-            H = a -> H_drift + sum(a .* H_drives)
-            G = a -> G_drift + sum(a .* G_drives)
-        end
+    return CompositeQuantumSystem{typeof(H), typeof(G)}(
+        H,
+        G,
+        H_drift,
+        H_drives,
+        T_max,
+        drive_bounds,
+        n_drives,
+        levels,
+        subsystem_levels,
+        subsystems
+    )
+end
 
-        return new{typeof(H), typeof(G)}(
-            H,
-            G,
-            n_drives,
-            levels,
-            params,
-            subsystem_levels,
-            subsystems
-        )
-    end
+"""
+    CompositeQuantumSystem(
+        H_drives::AbstractVector{<:AbstractMatrix},
+        subsystems::AbstractVector{<:QuantumSystem},
+        T_max::Float64,
+        drive_bounds::Vector
+    )
 
-    function CompositeQuantumSystem(
-        H_drives::AbstractVector{<:AbstractMatrix{T}},
-        subsystems::AbstractVector{<:QuantumSystem};
-        kwargs...
-    ) where T <: Number
-        @assert !isempty(H_drives) "At least one drive is required"
-        return CompositeQuantumSystem(
-            spzeros(T, size(H_drives[1])),
-            H_drives,
-            subsystems;
-            kwargs...
-        )
-    end
+Convenience constructor for a composite system with coupling drives but no coupling drift.
 
-    function CompositeQuantumSystem(
-        H_drift::AbstractMatrix{T},
-        subsystems::AbstractVector{<:QuantumSystem};
-        kwargs...
-    ) where T <: Number
-        return CompositeQuantumSystem(
-            H_drift, 
-            Matrix{T}[], 
-            subsystems;
-            kwargs...
-        )
-    end
+# Example
+```julia
+sys1 = QuantumSystem([PAULIS[:X]], 10.0, [1.0])
+sys2 = QuantumSystem([PAULIS[:Y]], 10.0, [1.0])
+g12 = 0.1 * kron(PAULIS[:X], PAULIS[:X])  # coupling drive
+csys = CompositeQuantumSystem([g12], [sys1, sys2], 10.0, [(-1.0, 1.0)])
+```
+"""
+function CompositeQuantumSystem(
+    H_drives::AbstractVector{<:AbstractMatrix{T}},
+    subsystems::AbstractVector{<:QuantumSystem},
+    T_max::Float64,
+    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}
+) where T <: Number
+    @assert !isempty(H_drives) "At least one drive is required"
+    return CompositeQuantumSystem(
+        spzeros(T, size(H_drives[1])),
+        H_drives,
+        subsystems,
+        T_max,
+        drive_bounds
+    )
+end
 
-    function CompositeQuantumSystem(
-        subsystems::AbstractVector{<:QuantumSystem};
-        kwargs...
+"""
+    CompositeQuantumSystem(
+        H_drift::AbstractMatrix,
+        subsystems::AbstractVector{<:QuantumSystem},
+        T_max::Float64,
+        drive_bounds::Vector
+    )
+
+Convenience constructor for a composite system with coupling drift but no coupling drives.
+
+# Example
+```julia
+sys1 = QuantumSystem([PAULIS[:X]], 10.0, [1.0])
+sys2 = QuantumSystem([PAULIS[:Y]], 10.0, [1.0])
+H_coupling = 0.1 * kron(PAULIS[:Z], PAULIS[:Z])  # coupling drift
+csys = CompositeQuantumSystem(H_coupling, [sys1, sys2], 10.0, Float64[])
+```
+"""
+function CompositeQuantumSystem(
+    H_drift::AbstractMatrix{T},
+    subsystems::AbstractVector{<:QuantumSystem},
+    T_max::Float64,
+    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}
+) where T <: Number
+    return CompositeQuantumSystem(
+        H_drift, 
+        Matrix{T}[], 
+        subsystems,
+        T_max,
+        drive_bounds
+    )
+end
+
+"""
+    CompositeQuantumSystem(
+        subsystems::AbstractVector{<:QuantumSystem},
+        T_max::Float64,
+        drive_bounds::Vector
+    )
+
+Convenience constructor for a composite system with no coupling terms (neither drift nor drives).
+
+Use this when you have independent subsystems that you want to represent in a single
+composite space, but without any direct coupling between them.
+
+# Example
+```julia
+sys1 = QuantumSystem([PAULIS[:X]], 10.0, [1.0])
+sys2 = QuantumSystem([PAULIS[:Y]], 10.0, [1.0])
+csys = CompositeQuantumSystem([sys1, sys2], 10.0, Float64[])
+```
+"""
+function CompositeQuantumSystem(
+    subsystems::AbstractVector{<:QuantumSystem},
+    T_max::Float64,
+    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}
+)
+    @assert !isempty(subsystems) "At least one subsystem is required"
+    T = eltype(get_drift(subsystems[1]))
+    levels = prod([sys.levels for sys ∈ subsystems])
+    return CompositeQuantumSystem(
+        spzeros(T, (levels, levels)), 
+        Matrix{T}[], 
+        subsystems,
+        T_max,
+        drive_bounds
     )
-        @assert !isempty(subsystems) "At least one subsystem is required"
-        T = eltype(get_drift(subsystems[1]))
-        levels = prod([sys.levels for sys ∈ subsystems])
-        return CompositeQuantumSystem(
-            spzeros(T, (levels, levels)), 
-            Matrix{T}[], 
-            subsystems; 
-            kwargs...
-        )
-    end
 end
 
 # ****************************************************************************** #
@@ -249,32 +378,33 @@ end
 
 @testitem "Composite system" begin
     subsystem_levels = [4, 2, 2]
-    sys1 = QuantumSystem(kron(PAULIS[:Z], PAULIS[:Z]), [kron(PAULIS[:X], PAULIS[:Y])])
-    sys2 = QuantumSystem([PAULIS[:Y], PAULIS[:Z]])
-    sys3 = QuantumSystem(zeros(2, 2))
+    sys1 = QuantumSystem(kron(PAULIS[:Z], PAULIS[:Z]), [kron(PAULIS[:X], PAULIS[:Y])], 1.0, [(-1.0, 1.0)])
+    sys2 = QuantumSystem([PAULIS[:Y], PAULIS[:Z]], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
+    sys3 = QuantumSystem(zeros(ComplexF64, 2, 2), 1.0)
     subsystems = [sys1, sys2, sys3]
     g12 = 0.1 * lift_operator([kron(PAULIS[:X], PAULIS[:X]), PAULIS[:X]], [1, 2], subsystem_levels)
     g23 = 0.2 * lift_operator([PAULIS[:Y], PAULIS[:Y]], [2, 3], subsystem_levels)
 
     # Construct composite system
-    csys = CompositeQuantumSystem(g12, [g23], [sys1, sys2, sys3])
+    csys = CompositeQuantumSystem(g12, [g23], [sys1, sys2, sys3], 1.0, [(-1.0, 1.0)])
     @test csys.levels == prod(subsystem_levels)
     @test csys.n_drives == 1 + sum([sys.n_drives for sys ∈ subsystems])
     @test csys.subsystems == subsystems
     @test csys.subsystem_levels == subsystem_levels
-    @test get_drift(csys) ≈ g12 + lift_operator(kron(PAULIS[:Z], PAULIS[:Z]), 1, subsystem_levels)end
+    @test get_drift(csys) ≈ g12 + lift_operator(kron(PAULIS[:Z], PAULIS[:Z]), 1, subsystem_levels)
+end
 
 @testitem "Composite system from drift" begin
     using LinearAlgebra
 
     subsystem_levels = [2, 2]
-    sys1 = QuantumSystem([PAULIS[:X], PAULIS[:Y]])
-    sys2 = QuantumSystem([PAULIS[:Y], PAULIS[:Z]])
+    sys1 = QuantumSystem([PAULIS[:X], PAULIS[:Y]], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
+    sys2 = QuantumSystem([PAULIS[:Y], PAULIS[:Z]], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
     subsystems = [sys1, sys2]
     g12 = 0.1 * kron(PAULIS[:X], PAULIS[:X])
 
     # Construct composite system from drift
-    csys = CompositeQuantumSystem(g12, [sys1, sys2])
+    csys = CompositeQuantumSystem(g12, [sys1, sys2], 1.0, Float64[])
     @test csys.levels == prod(subsystem_levels)
     @test csys.n_drives == sum([sys.n_drives for sys ∈ subsystems])
     @test csys.subsystems == subsystems
@@ -284,14 +414,14 @@ end
 
 @testitem "Composite system from drives" begin
     subsystem_levels = [2, 2, 2]
-    sys1 = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])
-    sys2 = QuantumSystem([PAULIS[:Y], PAULIS[:Z]])
-    sys3 = QuantumSystem(zeros(2, 2))
+    sys1 = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
+    sys2 = QuantumSystem([PAULIS[:Y], PAULIS[:Z]], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
+    sys3 = QuantumSystem(zeros(ComplexF64, 2, 2), 1.0)
     subsystems = [sys1, sys2, sys3]
     g12 = 0.1 * lift_operator([PAULIS[:X], PAULIS[:X]], [1, 2], subsystem_levels)
     g23 = 0.2 * lift_operator([PAULIS[:Y], PAULIS[:Y]], [2, 3], subsystem_levels)
 
-    csys = CompositeQuantumSystem([g12, g23], [sys1, sys2, sys3])
+    csys = CompositeQuantumSystem([g12, g23], [sys1, sys2, sys3], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
     @test csys.levels == prod(subsystem_levels)
     @test csys.n_drives == 2 + sum([sys.n_drives for sys ∈ subsystems])
     @test csys.subsystems == subsystems
diff --git a/src/quantum_systems/quantum_systems.jl b/src/quantum_systems/quantum_systems.jl
index 51a4381..0a5052d 100644
--- a/src/quantum_systems/quantum_systems.jl
+++ b/src/quantum_systems/quantum_systems.jl
@@ -1,3 +1,6 @@
+using ..Isomorphisms
+using SparseArrays: sparse, spzeros
+
 export QuantumSystem
 export OpenQuantumSystem
 export VariationalQuantumSystem
@@ -12,79 +15,168 @@ export VariationalQuantumSystem
 A struct for storing quantum dynamics.
 
 # Fields
-- `H::Function`: The Hamiltonian function, excluding dissipation: a -> H(a).
-- `G::Function`: The isomorphic generator function, including dissipation, a -> G(a).
-- `levels::Int`: The number of levels in the system.
-- `n_drives::Int`: The number of drives in the system.
-
+- `H::Function`: The Hamiltonian function: (u, t) -> H(u, t), where u is the control vector and t is time
+- `G::Function`: The isomorphic generator function: (u, t) -> G(u, t), including the Hamiltonian mapped to superoperator space
+- `H_drift::SparseMatrixCSC{ComplexF64, Int}`: The drift Hamiltonian (time-independent component)
+- `H_drives::Vector{SparseMatrixCSC{ComplexF64, Int}}`: The drive Hamiltonians (control-dependent components)
+- `T_max::Float64`: Maximum evolution time
+- `drive_bounds::Vector{Tuple{Float64, Float64}}`: Drive amplitude bounds for each control (lower, upper)
+- `n_drives::Int`: The number of control drives in the system
+- `levels::Int`: The number of levels (dimension) in the system
+
+See also [`OpenQuantumSystem`](@ref), [`VariationalQuantumSystem`](@ref).
 """
 struct QuantumSystem{F1<:Function, F2<:Function} <: AbstractQuantumSystem
     H::F1
     G::F2
+    H_drift::SparseMatrixCSC{ComplexF64, Int}
+    H_drives::Vector{SparseMatrixCSC{ComplexF64, Int}}
+    T_max::Float64
+    drive_bounds::Vector{Tuple{Float64, Float64}}
     n_drives::Int
     levels::Int
-    params::Dict{Symbol, Any}
+end
+
+"""
+    QuantumSystem(H::Function, T_max::Float64, drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}})
 
-    """
-        QuantumSystem(H_drift::Matrix{<:Number}, H_drives::Vector{Matrix{<:Number}}; kwargs...)
-        QuantumSystem(H_drift::Matrix{<:Number}; kwargs...)
-        QuantumSystem(H_drives::Vector{Matrix{<:Number}}; kwargs...)
-        QuantumSystem(H::Function, n_drives::Int; kwargs...)
+Construct a QuantumSystem from a Hamiltonian function.
 
-    Constructs a `QuantumSystem` object from the drift and drive Hamiltonian terms.
-    """
-    function QuantumSystem end
+# Arguments
+- `H::Function`: Hamiltonian function with signature (u, t) -> H(u, t) where u is the control vector and t is time
+- `T_max::Float64`: Maximum evolution time
+- `drive_bounds::Vector`: Drive amplitude bounds. Can be tuples `(lower, upper)` or scalars (interpreted as `(-bound, bound)`)
 
-    function QuantumSystem(
+# Example
+```julia
+# Define a time-dependent Hamiltonian
+H = (u, t) -> PAULIS[:Z] + u[1] * PAULIS[:X] + u[2] * PAULIS[:Y]
+sys = QuantumSystem(H, 10.0, [(-1.0, 1.0), (-1.0, 1.0)])
+```
+"""
+function QuantumSystem(
+    H::Function,
+    T_max::Float64,
+    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}
+)
+    drive_bounds = [
+        b isa Tuple ? b : (-b, b) for b in drive_bounds
+    ]
+
+    # Extract drift by evaluating with zero controls
+    H_drift = H(zeros(length(drive_bounds)), 0.0)
+    levels = size(H_drift, 1)
+
+    return QuantumSystem(
+        H,
+        (u, t) -> Isomorphisms.G(H(u, t)),
+        sparse(H_drift),
+        Vector{SparseMatrixCSC{ComplexF64, Int}}(),  # Empty drives vector for function-based systems
+        T_max,
+        drive_bounds,
+        length(drive_bounds),
+        levels
+    )
+end
+
+"""
+    QuantumSystem(
         H_drift::AbstractMatrix{<:Number},
-        H_drives::Vector{<:AbstractMatrix{<:Number}};
-        params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
+        H_drives::Vector{<:AbstractMatrix{<:Number}},
+        T_max::Float64,
+        drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}
     )
-        levels = size(H_drift, 1)
-        H_drift = sparse(H_drift)
-        G_drift = sparse(Isomorphisms.G(H_drift))
-
-        n_drives = length(H_drives)
-        H_drives = sparse.(H_drives)
-        G_drives = sparse.(Isomorphisms.G.(H_drives))
-
-        if n_drives == 0
-            H = a -> H_drift
-            G = a -> G_drift
-        else
-            H = a -> H_drift + sum(a .* H_drives)
-            G = a -> G_drift + sum(a .* G_drives)
-        end
-
-        return new{typeof(H), typeof(G)}(
-            H,
-            G,
-            n_drives,
-            levels,
-            params
-        )
-    end
 
-    function QuantumSystem(H_drives::Vector{<:AbstractMatrix{ℂ}}; kwargs...) where ℂ <: Number
-        @assert !isempty(H_drives) "At least one drive is required"
-        return QuantumSystem(spzeros(ℂ, size(H_drives[1])), H_drives; kwargs...)
+Construct a QuantumSystem from drift and drive Hamiltonian terms.
+
+# Arguments
+- `H_drift::AbstractMatrix`: The drift (time-independent) Hamiltonian
+- `H_drives::Vector{<:AbstractMatrix}`: Vector of drive Hamiltonians, one for each control
+- `T_max::Float64`: Maximum evolution time
+- `drive_bounds::Vector`: Drive amplitude bounds for each control. Can be tuples `(lower, upper)` or scalars
+
+The resulting Hamiltonian is: H(u, t) = H_drift + Σᵢ uᵢ * H_drives[i]
+
+# Example
+```julia
+sys = QuantumSystem(
+    PAULIS[:Z],                    # drift
+    [PAULIS[:X], PAULIS[:Y]],      # drives
+    10.0,                          # T_max
+    [(-1.0, 1.0), (-1.0, 1.0)]     # bounds
+)
+```
+"""
+function QuantumSystem(
+    H_drift::AbstractMatrix{<:Number},
+    H_drives::Vector{<:AbstractMatrix{<:Number}},
+    T_max::Float64,
+    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}
+)
+    drive_bounds = [
+        b isa Tuple ? b : (-b, b) for b in drive_bounds
+    ]
+
+    H_drift = sparse(H_drift)
+    G_drift = sparse(Isomorphisms.G(H_drift))
+
+    n_drives = length(H_drives)
+    H_drives = sparse.(H_drives)
+    G_drives = sparse.(Isomorphisms.G.(H_drives))
+
+    if n_drives == 0
+        H = (u, t) -> H_drift
+        G = (u, t) -> G_drift 
+    else
+        H = (u, t) -> H_drift + sum(u .* H_drives)
+        G = (u, t) -> G_drift + sum(u .* G_drives)
     end
 
-    QuantumSystem(H_drift::AbstractMatrix{ℂ}; kwargs...) where ℂ <: Number =
-        QuantumSystem(H_drift, Matrix{ℂ}[]; kwargs...)
+    levels = size(H_drift, 1)
+
+    return QuantumSystem(
+        H,
+        G,
+        H_drift,
+        H_drives,
+        T_max,
+        drive_bounds,
+        n_drives,
+        levels
+    )
+end
+
+# Convenience constructors
+"""
+    QuantumSystem(H_drives::Vector{<:AbstractMatrix}, T_max::Float64, drive_bounds::Vector)
 
-    function QuantumSystem(
-        H::F, 
-        n_drives::Int; 
-        params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
-    ) where F <: Function
-        G = a -> Isomorphisms.G(sparse(H(a)))
-        levels = size(H(zeros(n_drives)), 1)
-        return new{F, typeof(G)}(H, G, n_drives, levels, params)
-    end
+Convenience constructor for a system with no drift Hamiltonian (H_drift = 0).
+
+# Example
+```julia
+sys = QuantumSystem([PAULIS[:X], PAULIS[:Y]], 10.0, [1.0, 1.0])
+```
+"""
+function QuantumSystem(H_drives::Vector{<:AbstractMatrix{ℂ}}, T_max::Float64, drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}}) where ℂ <: Number
+    @assert !isempty(H_drives) "At least one drive is required"
+    return QuantumSystem(spzeros(ℂ, size(H_drives[1])), H_drives, T_max, drive_bounds)
+end
+
+"""
+    QuantumSystem(H_drift::AbstractMatrix, T_max::Float64)
 
+Convenience constructor for a system with only a drift Hamiltonian (no drives).
+
+# Example
+```julia
+sys = QuantumSystem(PAULIS[:Z], 10.0)
+```
+"""
+function QuantumSystem(H_drift::AbstractMatrix{ℂ}, T_max::Float64) where ℂ <: Number 
+    QuantumSystem(H_drift, Matrix{ℂ}[], T_max, Float64[])
 end
 
+
 # ----------------------------------------------------------------------------- #
 # OpenQuantumSystem
 # ----------------------------------------------------------------------------- #
@@ -94,120 +186,198 @@ end
 
 A struct for storing open quantum dynamics.
 
-# Additional fields
-- `dissipation_operators::Vector{AbstractMatrix}`: The dissipation operators.
+# Fields
+- `H::Function`: The Hamiltonian function: (u, t) -> H(u, t)
+- `𝒢::Function`: The Lindbladian generator function: u -> 𝒢(u)
+- `H_drift::SparseMatrixCSC{ComplexF64, Int}`: The drift Hamiltonian
+- `H_drives::Vector{SparseMatrixCSC{ComplexF64, Int}}`: The drive Hamiltonians
+- `T_max::Float64`: Maximum evolution time
+- `drive_bounds::Vector{Tuple{Float64, Float64}}`: Drive amplitude bounds
+- `n_drives::Int`: The number of control drives
+- `levels::Int`: The number of levels in the system
+- `dissipation_operators::Vector{SparseMatrixCSC{ComplexF64, Int}}`: The dissipation operators
+- `params::Dict{Symbol, Any}`: Additional parameters
 
 See also [`QuantumSystem`](@ref).
-
 """
 struct OpenQuantumSystem{F1<:Function, F2<:Function} <: AbstractQuantumSystem
     H::F1
     𝒢::F2
+    H_drift::SparseMatrixCSC{ComplexF64, Int}
+    H_drives::Vector{SparseMatrixCSC{ComplexF64, Int}}
+    T_max::Float64
+    drive_bounds::Vector{Tuple{Float64, Float64}}
     n_drives::Int
     levels::Int
-    dissipation_operators::Vector{Matrix{ComplexF64}}
+    dissipation_operators::Vector{SparseMatrixCSC{ComplexF64, Int}}
     params::Dict{Symbol, Any}
+end
 
-    """
+"""
     OpenQuantumSystem(
         H_drift::AbstractMatrix{<:Number},
-        H_drives::AbstractVector{<:AbstractMatrix{<:Number}}
-        dissipation_operators::AbstractVector{<:AbstractMatrix{<:Number}};
-        kwargs...
+        H_drives::AbstractVector{<:AbstractMatrix{<:Number}},
+        T_max::Float64,
+        drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}};
+        dissipation_operators::AbstractVector{<:AbstractMatrix{<:Number}}=Matrix{ComplexF64}[],
+        params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
     )
     OpenQuantumSystem(
-        H_drift::Matrix{<:Number}, 
-        H_drives::AbstractVector{Matrix{<:Number}};
+        H_drift::AbstractMatrix{<:Number}, 
+        T_max::Float64, 
+        drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}};
         dissipation_operators::AbstractVector{<:AbstractMatrix{<:Number}}=Matrix{ComplexF64}[],
-        kwargs...
+        params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
     )
-    OpenQuantumSystem(H_drift::Matrix{<:Number}; kwargs...)
-    OpenQuantumSystem(H_drives::Vector{Matrix{<:Number}}; kwargs...)
-    OpenQuantumSystem(H::Function, n_drives::Int; kwargs...)
-
-    Constructs an `OpenQuantumSystem` object from the drift and drive Hamiltonian terms and
-    dissipation operators.
-    """
-    function OpenQuantumSystem end
-
-    function OpenQuantumSystem(
-        H_drift::AbstractMatrix{<:Number},
-        H_drives::AbstractVector{<:AbstractMatrix{<:Number}};
+    OpenQuantumSystem(
+        H_drives::Vector{<:AbstractMatrix{<:Number}},
+        T_max::Float64, 
+        drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}};
         dissipation_operators::AbstractVector{<:AbstractMatrix{<:Number}}=Matrix{ComplexF64}[],
         params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
     )
-        levels = size(H_drift, 1)
-        H_drift = sparse(H_drift)
-        𝒢_drift = Isomorphisms.G(Isomorphisms.ad_vec(H_drift))
-
-        n_drives = length(H_drives)
-        H_drives = sparse.(H_drives)
-        𝒢_drives = Isomorphisms.G.(Isomorphisms.ad_vec.(H_drives))
-
-        if isempty(dissipation_operators)
-            𝒟 = spzeros(size(𝒢_drift))
-        else
-            𝒟 = sum(
-                Isomorphisms.iso_D(L) 
-                    for L ∈ sparse.(dissipation_operators)
-            )
-        end
-
-        if n_drives == 0
-            H = a -> H_drift
-            𝒢 = a -> 𝒢_drift + 𝒟
-        else
-            H = a -> H_drift + sum(a .* H_drives)
-            𝒢 = a -> 𝒢_drift + sum(a .* 𝒢_drives) + 𝒟
-        end
-
-        return new{typeof(H), typeof(𝒢)}(
-            H,
-            𝒢,
-            n_drives,
-            levels,
-            dissipation_operators,
-            params
-        )
-    end
-
-    function OpenQuantumSystem(
-        H_drift::AbstractMatrix{<:Number},
-        H_drives::AbstractVector{<:AbstractMatrix{<:Number}},
-        dissipation_operators::AbstractVector{<:AbstractMatrix{<:Number}};
+    OpenQuantumSystem(
+        H::Function, 
+        T_max::Float64,
+        drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}};
+        dissipation_operators::Vector{<:AbstractMatrix{<:Number}}=Matrix{ComplexF64}[],
         params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
     )
-        return OpenQuantumSystem(
-            H_drift, H_drives;
-            dissipation_operators=dissipation_operators,
-            params=params
-        )
+    OpenQuantumSystem(
+        system::QuantumSystem; 
+        dissipation_operators::Vector{<:AbstractMatrix{<:Number}}=Matrix{ComplexF64}[],
+        params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
+    )
+
+Constructs an OpenQuantumSystem object from the drift and drive Hamiltonian terms and
+dissipation operators. All constructors require T_max (maximum time) and drive_bounds
+(control bounds for each drive) to be explicitly specified.
+"""
+function OpenQuantumSystem end
+
+# Main constructor from Hamiltonian components 
+function OpenQuantumSystem(
+    H_drift::AbstractMatrix{<:Number},
+    H_drives::Vector{<:AbstractMatrix{<:Number}},
+    T_max::Float64,
+    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}};
+    dissipation_operators::Vector{<:AbstractMatrix{<:Number}}=Matrix{ComplexF64}[],
+    params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
+)
+    drive_bounds = [
+        b isa Tuple ? b : (-b, b) for b in drive_bounds
+    ]
+
+    H_drift_sparse = sparse(H_drift)
+    𝒢_drift = Isomorphisms.G(Isomorphisms.ad_vec(H_drift_sparse))
+
+    n_drives = length(H_drives)
+    H_drives_sparse = sparse.(H_drives)
+    𝒢_drives = [Isomorphisms.G(Isomorphisms.ad_vec(H_drive)) for H_drive in H_drives_sparse]
+    
+    # Build dissipator
+    if isempty(dissipation_operators)
+        𝒟 = spzeros(size(𝒢_drift))
+    else
+        𝒟 = sum(Isomorphisms.iso_D(sparse(L)) for L in dissipation_operators)
     end
 
-    function OpenQuantumSystem(
-        H_drives::AbstractVector{<:AbstractMatrix{ℂ}}; kwargs...
-    ) where ℂ <: Number
-        @assert !isempty(H_drives) "At least one drive is required"
-        return OpenQuantumSystem(spzeros(ℂ, size(H_drives[1])), H_drives; kwargs...)
+    if n_drives == 0
+        H = (u, t) -> H_drift_sparse
+        𝒢 = u -> 𝒢_drift + 𝒟
+    else
+        H = (u, t) -> H_drift_sparse + sum(u .* H_drives_sparse)
+        𝒢 = u -> 𝒢_drift + sum(u .* 𝒢_drives) + 𝒟
     end
 
-    OpenQuantumSystem(H_drift::AbstractMatrix{T}; kwargs...) where T <: Number =
-        OpenQuantumSystem(H_drift, Matrix{T}[]; kwargs...)
+    levels = size(H_drift, 1)
+
+    return OpenQuantumSystem(
+        H,
+        𝒢,
+        H_drift_sparse,
+        H_drives_sparse,
+        T_max,
+        drive_bounds,
+        n_drives,
+        levels,
+        sparse.(dissipation_operators),
+        params
+    )
+end
 
-    function OpenQuantumSystem(
-        H::F, n_drives::Int;
-        dissipation_operators::AbstractVector{<:AbstractMatrix{ℂ}}=Matrix{ComplexF64}[],
-        params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
-    ) where {F <: Function, ℂ <: Number}
-        G = a -> Isomorphisms.G(Isomorphisms.ad_vec(sparse(H(a))))
-        levels = size(H(zeros(n_drives)), 1)
-        return new{F, typeof(G)}(H, G, n_drives, levels, dissipation_operators, params)
+# Convenience constructors
+function OpenQuantumSystem(
+    H_drives::Vector{<:AbstractMatrix{ℂ}}, 
+    T_max::Float64, 
+    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}};
+    dissipation_operators::Vector{<:AbstractMatrix{<:Number}}=Matrix{ComplexF64}[],
+    params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
+) where ℂ <: Number
+    @assert !isempty(H_drives) "At least one drive is required"
+    return OpenQuantumSystem(spzeros(ℂ, size(H_drives[1])), H_drives, T_max, drive_bounds;
+                            dissipation_operators=dissipation_operators, params=params)
+end
+
+function OpenQuantumSystem(
+    H_drift::AbstractMatrix{ℂ}, 
+    T_max::Float64; 
+    dissipation_operators::Vector{<:AbstractMatrix{<:Number}}=Matrix{ComplexF64}[],
+    params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
+) where ℂ <: Number 
+    return OpenQuantumSystem(H_drift, Matrix{ℂ}[], T_max, Float64[];
+                            dissipation_operators=dissipation_operators, params=params)
+end
+
+function OpenQuantumSystem(
+    H::F, 
+    T_max::Float64,
+    drive_bounds::Vector{<:Union{Tuple{Float64, Float64}, Float64}};
+    dissipation_operators::Vector{<:AbstractMatrix{ℂ}}=Matrix{ComplexF64}[],
+    params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
+) where {F <: Function, ℂ <: Number}
+    
+    drive_bounds = [
+        b isa Tuple ? b : (-b, b) for b in drive_bounds
+    ]
+
+    n_drives = length(drive_bounds)
+    
+    # Extract drift by evaluating with zero controls
+    H_drift = H(zeros(n_drives), 0.0)
+    levels = size(H_drift, 1)
+    
+    # Build dissipator
+    if isempty(dissipation_operators)
+        𝒟 = spzeros(ComplexF64, levels^2, levels^2)
+    else
+        𝒟 = sum(Isomorphisms.iso_D(sparse(L)) for L in dissipation_operators)
     end
 
-    OpenQuantumSystem(system::QuantumSystem; kwargs...) = OpenQuantumSystem(
-        system.H, system.n_drives; kwargs...
+    return OpenQuantumSystem(
+        H,
+        u -> Isomorphisms.G(Isomorphisms.ad_vec(sparse(H(u, 0.0)))) + 𝒟,
+        sparse(H_drift),
+        Vector{SparseMatrixCSC{ComplexF64, Int}}(),  # Empty drives vector for function-based systems
+        T_max,
+        drive_bounds,
+        n_drives,
+        levels,
+        sparse.(dissipation_operators),
+        params
     )
+end
 
+function OpenQuantumSystem(
+    system::QuantumSystem; 
+    dissipation_operators::Vector{<:AbstractMatrix{<:Number}}=Matrix{ComplexF64}[],
+    params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
+)
+    return OpenQuantumSystem(
+        system.H_drift, system.H_drives, system.T_max, system.drive_bounds;
+        dissipation_operators=dissipation_operators,
+        params=params
+    )
 end
 
 # ----------------------------------------------------------------------------- #
@@ -219,13 +389,20 @@ end
 """
     VariationalQuantumSystem <: AbstractQuantumSystem
 
-A struct for storing variational quantum dynamics.
+A struct for storing variational quantum dynamics, used for sensitivity and robustness analysis.
 
-# Additional fields
-- `G_vars::AbstractVector{<:Function}`: Variational generator functions
-
-See also [`QuantumSystem`](@ref).
+Variational systems allow exploring how the dynamics change under perturbations to the Hamiltonian.
+The variational operators represent directions of uncertainty or perturbation in the system.
 
+# Fields
+- `H::Function`: The Hamiltonian function: (u, t) -> H(u, t)
+- `G::Function`: The isomorphic generator function: (u, t) -> G(u, t)
+- `G_vars::AbstractVector{<:Function}`: Variational generator functions, one for each perturbation direction
+- `n_drives::Int`: The number of control drives in the system
+- `levels::Int`: The number of levels (dimension) in the system
+- `params::Dict{Symbol, Any}`: Additional parameters for the system
+
+See also [`QuantumSystem`](@ref), [`OpenQuantumSystem`](@ref).
 """
 struct VariationalQuantumSystem{F1<:Function, F2<:Function, F⃗3<:AbstractVector{<:Function}} <: AbstractQuantumSystem
     H::F1 
@@ -234,172 +411,199 @@ struct VariationalQuantumSystem{F1<:Function, F2<:Function, F⃗3<:AbstractVecto
     n_drives::Int 
     levels::Int 
     params::Dict{Symbol, Any}
+end
 
-    """
-        VariationalQuantumSystem(
-            H_drift::AbstractMatrix{<:Number},
-            H_drives::AbstractVector{<:AbstractMatrix{<:Number}},
-            H_vars::AbstractVector{<:AbstractMatrix{<:Number}};
-            kwargs...
-        )
-        VariationalQuantumSystem(
-            H_drives::AbstractVector{<:AbstractMatrix{<:Number}},
-            H_vars::AbstractVector{<:AbstractMatrix{<:Number}};
-            kwargs...
-        )
-        VariationalQuantumSystem(
-            H::F1,
-            H_vars::F⃗2,
-            n_drives::Int;
-            kwargs...
-        )
-
-    Constructs a `VariationalQuantumSystem` object from the drift and drive Hamiltonian 
-    terms and variational Hamiltonians, for sensitivity and robustness analysis.
-
-    """
-    function VariationalQuantumSystem end
-
-    function VariationalQuantumSystem(
-        H_drift::AbstractMatrix{<:Number},
-        H_drives::AbstractVector{<:AbstractMatrix{<:Number}},
-        H_vars::AbstractVector{<:AbstractMatrix{<:Number}};
-        params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
+"""
+    VariationalQuantumSystem(
+        H_drift::AbstractMatrix,
+        H_drives::AbstractVector{<:AbstractMatrix},
+        H_vars::AbstractVector{<:AbstractMatrix};
+        params::Dict{Symbol,Any}=Dict{Symbol,Any}()
     )
-        @assert !isempty(H_vars) "At least one variational operator is required"
 
-        levels = size(H_drift, 1)
-        H_drift = sparse(H_drift)
-        G_drift = sparse(Isomorphisms.G(H_drift))
+Construct a VariationalQuantumSystem from drift, drive, and variational Hamiltonian terms.
+
+# Arguments
+- `H_drift::AbstractMatrix`: The drift (time-independent) Hamiltonian
+- `H_drives::AbstractVector{<:AbstractMatrix}`: Vector of drive Hamiltonians for control
+- `H_vars::AbstractVector{<:AbstractMatrix}`: Vector of variational Hamiltonians representing perturbation directions
+- `params::Dict{Symbol,Any}`: Optional additional parameters
+
+The variational operators allow sensitivity analysis by exploring how dynamics change
+under perturbations: H_perturbed = H + Σᵢ εᵢ * H_vars[i]
+
+# Example
+```julia
+varsys = VariationalQuantumSystem(
+    PAULIS[:Z],                    # drift
+    [PAULIS[:X], PAULIS[:Y]],      # drives
+    [PAULIS[:X]]                   # variational perturbations
+)
+```
+"""
+function VariationalQuantumSystem end
 
-        n_drives = length(H_drives)
-        H_drives = sparse.(H_drives)
-        G_drives = sparse.(Isomorphisms.G.(H_drives))
+function VariationalQuantumSystem(
+    H_drift::AbstractMatrix{<:Number},
+    H_drives::AbstractVector{<:AbstractMatrix{<:Number}},
+    H_vars::AbstractVector{<:AbstractMatrix{<:Number}};
+    params::Dict{Symbol,Any}=Dict{Symbol,Any}()
+)
+    @assert !isempty(H_vars) "At least one variational operator is required"
 
-        G_vars = [a -> Isomorphisms.G(sparse(H)) for H in H_vars]
+    levels = size(H_drift, 1)
+    H_drift = sparse(H_drift)
+    G_drift = sparse(Isomorphisms.G(H_drift))
 
-        if n_drives == 0
-            H = a -> H_drift
-            G = a -> G_drift
-        else
-            H = a -> H_drift + sum(a .* H_drives)
-            G = a -> G_drift + sum(a .* G_drives)
-        end
+    n_drives = length(H_drives)
+    H_drives = sparse.(H_drives)
+    G_drives = sparse.(Isomorphisms.G.(H_drives))
 
-        return new{typeof(H), typeof(G), typeof(G_vars)}(
-            H, G, G_vars, n_drives, levels, params
-        )
-    end
+    G_vars = [a -> Isomorphisms.G(sparse(H)) for H in H_vars]
 
-    function VariationalQuantumSystem(
-        H_drives::AbstractVector{<:AbstractMatrix{ℂ}},
-        H_vars::AbstractVector{<:AbstractMatrix{<:Number}};
-        kwargs...
-    ) where ℂ <: Number
-        @assert !isempty(H_drives) "At least one drive is required"
-        @assert !isempty(H_vars) "At least one variational operator is required"
-        return VariationalQuantumSystem(
-            spzeros(ℂ, size(H_drives[1])), 
-            H_drives, 
-            H_vars; 
-            kwargs...
-        )
+    if n_drives == 0
+        H = a -> H_drift
+        G = a -> G_drift
+    else
+        H = a -> H_drift + sum(a .* H_drives)
+        G = a -> G_drift + sum(a .* G_drives)
     end
 
-    function VariationalQuantumSystem(
-        H::F1,
-        H_vars::F⃗2,
-        n_drives::Int; 
-        params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
-    ) where {F1 <: Function, F⃗2 <: AbstractVector{<:Function}}
-        @assert !isempty(H_vars) "At least one variational operator is required"
-        G = a -> Isomorphisms.G(sparse(H(a)))
-        G_vars = Function[a -> Isomorphisms.G(sparse(H_v(a))) for H_v in H_vars]
-        levels = size(H(zeros(n_drives)), 1)
-        return new{F1, typeof(G), F⃗2}(H, G, G_vars, n_drives, levels, params)
-    end
+    return VariationalQuantumSystem(
+        H, G, G_vars, n_drives, levels, params
+    )
+end
+
+function VariationalQuantumSystem(
+    H_drives::AbstractVector{<:AbstractMatrix{ℂ}},
+    H_vars::AbstractVector{<:AbstractMatrix{<:Number}};
+    params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
+) where ℂ <: Number
+    @assert !isempty(H_drives) "At least one drive is required"
+    @assert !isempty(H_vars) "At least one variational operator is required"
+    return VariationalQuantumSystem(
+        spzeros(ℂ, size(H_drives[1])), 
+        H_drives, 
+        H_vars; 
+        params=params
+    )
+end
+
+function VariationalQuantumSystem(
+    H::F1,
+    H_vars::F⃗2,
+    n_drives::Int; 
+    params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
+) where {F1 <: Function, F⃗2 <: AbstractVector{<:Function}}
+    @assert !isempty(H_vars) "At least one variational operator is required"
+    G = a -> Isomorphisms.G(sparse(H(a)))
+    G_vars = Function[a -> Isomorphisms.G(sparse(H_v(a))) for H_v in H_vars]
+    levels = size(H(zeros(n_drives)), 1)
+    return VariationalQuantumSystem(H, G, G_vars, n_drives, levels, params)
 end
 
 # ******************************************************************************* #
 
 @testitem "System creation" begin
-    H_drift = PAULIS[:Z]
-    H_drives = [PAULIS[:X], PAULIS[:Y]]
+    using PiccoloQuantumObjects: PAULIS, QuantumSystem, get_drift, get_drives
+    using SparseArrays: sparse
+    
+    H_drift = PAULIS.Z
+    H_drives = [PAULIS.X, PAULIS.Y]
     n_drives = length(H_drives)
-
-    system = QuantumSystem(H_drift, H_drives)
+    T_max = 1.0
+    u_bounds = ones(n_drives)
+    
+    system = QuantumSystem(H_drift, H_drives, T_max, u_bounds)
     @test system isa QuantumSystem
     @test get_drift(system) == H_drift
     @test get_drives(system) == H_drives
 
     # repeat with a bigger system
-    H_drift = kron(PAULIS[:Z], PAULIS[:Z])
-    H_drives = [kron(PAULIS[:X], PAULIS[:I]), kron(PAULIS[:I], PAULIS[:X]),
-                kron(PAULIS[:Y], PAULIS[:I]), kron(PAULIS[:I], PAULIS[:Y])]
+    H_drift = kron(PAULIS.Z, PAULIS.Z)
+    H_drives = [kron(PAULIS.X, PAULIS.I), kron(PAULIS.I, PAULIS.X),
+                kron(PAULIS.Y, PAULIS.I), kron(PAULIS.I, PAULIS.Y)]
     n_drives = length(H_drives)
+    u_bounds = ones(n_drives)
 
-    system = QuantumSystem(H_drift, H_drives)
-    @test system isa AbstractQuantumSystem
+    system = QuantumSystem(H_drift, H_drives, T_max, u_bounds)
     @test system isa QuantumSystem
     @test get_drift(system) == H_drift
     @test get_drives(system) == H_drives
 end
 
-@testitem "System creation with params" begin
-    system = QuantumSystem(PAULIS[:Z], [PAULIS[:X]], params=Dict(:a => 1))
-    @test system.params[:a] == 1
-
-    open_system = OpenQuantumSystem(PAULIS[:Z], [PAULIS[:X]], params=Dict(:a => 1))
-    @test open_system.params[:a] == 1
-end
-
 @testitem "No drift system creation" begin
-    H_drift = zeros(2, 2)
-    H_drives = [PAULIS[:X], PAULIS[:Y]]
+    using PiccoloQuantumObjects: PAULIS, QuantumSystem, get_drift, get_drives
+    using SparseArrays: spzeros
+    
+    H_drift = zeros(ComplexF64, 2, 2)
+    H_drives = [PAULIS.X, PAULIS.Y]
+    T_max = 1.0
+    u_bounds = [1.0, 1.0]
 
-    sys1 = QuantumSystem(H_drift, H_drives)
-    sys2 = QuantumSystem(H_drives)
+    sys1 = QuantumSystem(H_drift, H_drives, T_max, u_bounds)
+    sys2 = QuantumSystem(H_drives, T_max, u_bounds)
 
     @test get_drift(sys1) == get_drift(sys2) == H_drift
     @test get_drives(sys1) == get_drives(sys2) == H_drives
 end
 
 @testitem "No drive system creation" begin
-    H_drift = PAULIS[:Z]
+    using PiccoloQuantumObjects: PAULIS, QuantumSystem, get_drift, get_drives
+    
+    H_drift = PAULIS.Z
     H_drives = Matrix{ComplexF64}[]
+    T_max = 1.0
+    u_bounds = Float64[]
 
-    sys1 = QuantumSystem(H_drift, H_drives)
-    sys2 = QuantumSystem(H_drift)
+    sys1 = QuantumSystem(H_drift, H_drives, T_max, u_bounds)
+    sys2 = QuantumSystem(H_drift, T_max)
 
     @test get_drift(sys1) == get_drift(sys2) == H_drift
     @test get_drives(sys1) == get_drives(sys2) == H_drives
 end
 
 @testitem "System creation with Hamiltonian function" begin
+    using PiccoloQuantumObjects: PAULIS, QuantumSystem, get_drift, get_drives
+
+    # test one drive
+
+    H_drift = PAULIS.Z
+    H_drives = [PAULIS.X]
+    
     system = QuantumSystem(
-        a -> PAULIS[:Z] + a[1] * PAULIS[:X], 1
+        (a, t) -> H_drift + sum(a .* H_drives), 
+        1.0, 
+        [1.0]
     )
     @test system isa QuantumSystem
-    @test get_drift(system) == PAULIS[:Z]
-    @test get_drives(system) == [PAULIS[:X]]
+    @test get_drift(system) == H_drift 
+    @test get_drives(system) == H_drives
+
+    # test no drift + three drives
 
-    # test three drives
+    H_drives = [PAULIS.X, PAULIS.Y, PAULIS.Z]
     system = QuantumSystem(
-        a -> a[1] * PAULIS[:X] + a[2] * PAULIS[:Y] + a[3] * PAULIS[:Z], 3
+        (a, t) -> sum(a .* H_drives),
+        1.0, 
+        [1.0, 1.0, 1.0]
     )
     @test system isa QuantumSystem
     @test get_drift(system) == zeros(2, 2)
-    @test get_drives(system) == [PAULIS[:X], PAULIS[:Y], PAULIS[:Z]]
+    @test get_drives(system) == H_drives 
 end
 
 @testitem "Open system creation" begin
-    H_drift = PAULIS[:Z]
+    using PiccoloQuantumObjects: PAULIS, OpenQuantumSystem, get_drift, get_drives, Isomorphisms
+    
+    H_drift = PAULIS.Z
     # don't want drives == levels
-    H_drives = [PAULIS[:X]]
-    dissipation_operators = [PAULIS[:Z], PAULIS[:X]]
+    H_drives = [PAULIS.X]
+    dissipation_operators = [PAULIS.Z, PAULIS.X]
+    T_max = 1.0
+    drive_bounds = [1.0]
 
-    system = OpenQuantumSystem(H_drift, H_drives, dissipation_operators)
-    @test system isa AbstractQuantumSystem
+    system = OpenQuantumSystem(H_drift, H_drives, T_max, drive_bounds, dissipation_operators=dissipation_operators)
     @test system isa OpenQuantumSystem
     @test get_drift(system) == H_drift
     @test get_drives(system) == H_drives
@@ -411,13 +615,17 @@ end
 end
 
 @testitem "Open system alternate constructors" begin
-    H_drift = PAULIS[:Z]
+    using PiccoloQuantumObjects: PAULIS, OpenQuantumSystem, get_drift, get_drives
+    
+    H_drift = PAULIS.Z
     # don't want drives == levels
-    H_drives = [PAULIS[:X]]
-    dissipation_operators = [PAULIS[:Z], PAULIS[:X]]
+    H_drives = [PAULIS.X]
+    dissipation_operators = [PAULIS.Z, PAULIS.X]
+    T_max = 1.0
+    drive_bounds = [1.0]
 
     system = OpenQuantumSystem(
-        H_drift, H_drives, dissipation_operators=dissipation_operators
+        H_drift, H_drives, T_max, drive_bounds, dissipation_operators=dissipation_operators
     )
     @test system isa OpenQuantumSystem
     @test get_drift(system) == H_drift
@@ -425,7 +633,7 @@ end
     @test system.dissipation_operators == dissipation_operators
 
     # no drift
-    system = OpenQuantumSystem(H_drives, dissipation_operators=dissipation_operators)
+    system = OpenQuantumSystem(H_drives, T_max, drive_bounds, dissipation_operators=dissipation_operators)
     @test system isa OpenQuantumSystem
     @test get_drift(system) == zeros(size(H_drift))
     @test get_drives(system) == H_drives
@@ -433,7 +641,7 @@ end
 
     # no drives
     system = OpenQuantumSystem(
-        H_drift, dissipation_operators=dissipation_operators
+        H_drift, T_max, dissipation_operators=dissipation_operators
     )
     @test system isa OpenQuantumSystem
     @test system isa OpenQuantumSystem
@@ -442,8 +650,16 @@ end
     @test system.dissipation_operators == dissipation_operators
 
     # function
-    H = a -> PAULIS[:Z] + a[1] * PAULIS[:X]
-    system = OpenQuantumSystem(H, 1, dissipation_operators=dissipation_operators)
+    H = (u, t) -> PAULIS.Z + u[1] * PAULIS.X
+    system = OpenQuantumSystem(H, T_max, drive_bounds, dissipation_operators=dissipation_operators)
+    @test system isa OpenQuantumSystem
+    @test get_drift(system) == H_drift
+    @test get_drives(system) == H_drives
+    @test system.dissipation_operators == dissipation_operators
+
+    # from QuantumSystem
+    qsys = QuantumSystem(H_drift, H_drives, T_max, drive_bounds)
+    system = OpenQuantumSystem(qsys, dissipation_operators=dissipation_operators)
     @test system isa OpenQuantumSystem
     @test get_drift(system) == H_drift
     @test get_drives(system) == H_drives
@@ -452,6 +668,9 @@ end
 end
 
 @testitem "Variational system creation" begin
+    using PiccoloQuantumObjects: PAULIS, VariationalQuantumSystem, Isomorphisms
+    using LinearAlgebra: I
+    
     # default
     varsys1 = VariationalQuantumSystem(
         0.0 * PAULIS.Z,
@@ -470,12 +689,12 @@ end
     G_Y = Isomorphisms.G(PAULIS.Y)
     G = a[1] * G_X + a[2] * G_Y
     for varsys in [varsys1, varsys2]
-        @assert varsys isa VariationalQuantumSystem
-        @assert varsys.n_drives == 2
-        @assert length(varsys.G_vars) == 2
-        @assert varsys.G(a) ≈ G
-        @assert varsys.G_vars[1](a) ≈ G_X
-        @assert varsys.G_vars[2](a) ≈ G_Y
+        @test varsys isa VariationalQuantumSystem
+        @test varsys.n_drives == 2
+        @test length(varsys.G_vars) == 2
+        @test varsys.G(a) ≈ G
+        @test varsys.G_vars[1](a) ≈ G_X
+        @test varsys.G_vars[2](a) ≈ G_Y
     end
 
     # single sensitivity
@@ -483,11 +702,11 @@ end
         [PAULIS.X, PAULIS.Y],
         [PAULIS.X] 
     )
-    @assert varsys isa VariationalQuantumSystem
-    @assert varsys.n_drives == 2
-    @assert length(varsys.G_vars) == 1
-    @assert varsys.G(a) ≈ G
-    @assert varsys.G_vars[1](a) ≈ G_X
+    @test varsys isa VariationalQuantumSystem
+    @test varsys.n_drives == 2
+    @test length(varsys.G_vars) == 1
+    @test varsys.G(a) ≈ G
+    @test varsys.G_vars[1](a) ≈ G_X
 
     # functional sensitivity
     varsys = VariationalQuantumSystem(
@@ -495,10 +714,10 @@ end
         [a -> a[1] * PAULIS.X, a -> PAULIS.Y],
         2
     )
-    @assert varsys isa VariationalQuantumSystem
-    @assert varsys.n_drives == 2
-    @assert length(varsys.G_vars) == 2
-    @assert varsys.G(a) ≈ G
-    @assert varsys.G_vars[1](a) ≈ a[1] * G_X
-    @assert varsys.G_vars[2](a) ≈ G_Y
+    @test varsys isa VariationalQuantumSystem
+    @test varsys.n_drives == 2
+    @test length(varsys.G_vars) == 2
+    @test varsys.G(a) ≈ G
+    @test varsys.G_vars[1](a) ≈ a[1] * G_X
+    @test varsys.G_vars[2](a) ≈ G_Y
 end
\ No newline at end of file
diff --git a/src/quantum_systems/time_dependent_quantum_systems.jl b/src/quantum_systems/time_dependent_quantum_systems.jl
deleted file mode 100644
index ff8811f..0000000
--- a/src/quantum_systems/time_dependent_quantum_systems.jl
+++ /dev/null
@@ -1,200 +0,0 @@
-export TimeDependentQuantumSystem
-
-# ----------------------------------------------------------------------------- #
-# TimeDependentQuantumSystem
-# ----------------------------------------------------------------------------- #
-
-# TODO: Open System
-
-"""
-    TimeDependentQuantumSystem <: AbstractQuantumSystem
-
-A struct for storing time-dependent quantum dynamics.
-
-# Fields
-- `H::Function`: The Hamiltonian function with time: (a, t) -> H(a, t).
-- `G::Function`: The isomorphic generator function with time, (a, t) -> G(a, t).
-- `n_drives::Int`: The number of drives in the system.
-- `levels::Int`: The number of levels in the system.
-- `params::Dict{Symbol, Any}`: A dictionary of parameters.
-
-"""
-struct TimeDependentQuantumSystem{F1<:Function, F2<:Function} <: AbstractQuantumSystem
-    H::F1
-    G::F2
-    n_drives::Int
-    levels::Int 
-    params::Dict{Symbol, Any}
-
-    """
-        TimeDependentQuantumSystem(
-            H::AbstractMatrix{<:Number}, 
-            H_drives::Vector{<:AbstractMatrix{<:Number}},
-            carriers::Vector{<:Real}, 
-            phases::Vector{<:Real}; 
-            kwargs...
-        )
-        TimeDependentQuantumSystem(
-            H::AbstractMatrix{<:Number}, 
-            H_drives::Vector{<:AbstractMatrix{<:Number}},
-            kwargs...
-        )
-        TimeDependentQuantumSystem(
-            H::Function,
-            n_drives::Int;
-            params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
-        )
-
-    Constructs a time-dependent quantum system with the given Hamiltonian, drives,
-    carriers, and phases.
-    """
-    function TimeDependentQuantumSystem end
-
-    function TimeDependentQuantumSystem(
-        H_drift::AbstractMatrix{<:Number},
-        H_drives::Vector{<:AbstractMatrix{<:Number}};
-        carriers::AbstractVector{<:Real}=zeros(length(H_drives)),
-        phases::AbstractVector{<:Real}=zeros(length(H_drives)),
-        params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
-    )
-        levels = size(H_drift, 1)
-        H_drift = sparse(H_drift)
-        G_drift = sparse(Isomorphisms.G(H_drift))
-
-        n_drives = length(H_drives)
-        H_drives = sparse.(H_drives)
-        G_drives = sparse.(Isomorphisms.G.(H_drives))
-
-        if n_drives == 0
-            H = (a, t) -> H_drift
-            G = (a, t) -> G_drift
-        else
-            H = (a, t) -> H_drift + sum(a[i] * cos(carriers[i] * t + phases[i]) .* H_drives[i] for i in eachindex(H_drives))
-            G = (a, t) -> G_drift + sum(a[i] * cos(carriers[i] * t + phases[i]) .* G_drives[i] for i in eachindex(G_drives))
-        end
-
-        # save carries and phases
-        params[:carriers] = carriers
-        params[:phases] = phases
-
-        return new{typeof(H), typeof(G)}(
-            H,
-            G,
-            n_drives,
-            levels,
-            params
-        )
-    end
-
-    function TimeDependentQuantumSystem(
-        H_drives::Vector{<:AbstractMatrix{ℂ}}; 
-        kwargs...
-    ) where ℂ <: Number
-        @assert !isempty(H_drives) "At least one drive is required."
-        return TimeDependentQuantumSystem(
-            spzeros(ℂ, size(H_drives[1])), 
-            H_drives; 
-            kwargs...
-        )
-    end
-
-    TimeDependentQuantumSystem(H_drift::AbstractMatrix{ℂ}; kwargs...) where ℂ <: Number = 
-        TimeDependentQuantumSystem(H_drift, Matrix{ℂ}[]; kwargs...)
-
-    function TimeDependentQuantumSystem(
-        H::F,
-        n_drives::Int;
-        params::Dict{Symbol, <:Any}=Dict{Symbol, Any}()
-    ) where F <: Function
-        G = (a, t) -> Isomorphisms.G(sparse(H(a, t)))
-        levels = size(H(zeros(n_drives), 0.0), 1)
-        return new{F, typeof(G)}(H, G, n_drives, levels, params)
-    end
-end
-
-get_drift(sys::TimeDependentQuantumSystem) = t -> sys.H(zeros(sys.n_drives), t)
-
-function get_drives(sys::TimeDependentQuantumSystem)
-    H_drift = get_drift(sys)
-    # Basis vectors for controls will extract drive operators
-    return [t -> sys.H(I[1:sys.n_drives, i], t) - H_drift(t) for i ∈ 1:sys.n_drives]
-end
-
-# ****************************************************************************** #
-
-@testitem "System creation" begin
-    H_drift = kron(PAULIS.Z, PAULIS.I)
-    H_drives = [kron(PAULIS.X, PAULIS.I), kron(PAULIS.Y, PAULIS.I)]
-    n_drives = length(H_drives)
-    carriers = [1.0, 2.0]
-    phases = [0.1, 0.2]
-
-    sys = TimeDependentQuantumSystem(H_drift, H_drives, carriers=carriers, phases=phases)
-
-    @test sys isa TimeDependentQuantumSystem
-
-    # Params
-    @test sys.params[:carriers] == carriers
-    @test sys.params[:phases] == phases
-
-    # Time-dependent Hamiltonians
-    @test get_drift(sys)(0.0) == H_drift
-    @test [H(0.0) for H in get_drives(sys)] == H_drives .* cos.(phases)
-
-    t = 2.0
-    @test get_drift(sys)(t) == H_drift
-    @test [H(t) for H in get_drives(sys)] == H_drives .* cos.(carriers .* t .+ phases)
-end
-
-@testitem "No drift system creation" begin
-    H_drift = zeros(2, 2)
-    H_drives = [PAULIS.X, PAULIS.Y]
-    sys1 = TimeDependentQuantumSystem(H_drift, H_drives)
-    sys2 = TimeDependentQuantumSystem(H_drives)
-
-    @test get_drift(sys1)(0.0) == get_drift(sys2)(0.0) == H_drift
-    @test [H(0.0) for H in get_drives(sys1)] == [H(0.0) for H in get_drives(sys2)] == H_drives
-end
-
-@testitem "No drive system creation" begin
-    H_drift = PAULIS.Z
-    H_drives = Matrix{ComplexF64}[]
-
-    sys1 = TimeDependentQuantumSystem(H_drift, H_drives)
-    sys2 = TimeDependentQuantumSystem(H_drift)
-
-    @test get_drift(sys1)(0.0) == get_drift(sys2)(0.0) == H_drift
-    @test [H(0.0) for H in get_drives(sys1)] == [H(0.0) for H in get_drives(sys2)] == H_drives
-end
-
-@testitem "System creation with Hamiltonian function" begin
-    n_drives = 1
-    H_drift = t -> PAULIS.Z * cos(t)
-    H_drive = t -> PAULIS.X * cos(t)
-    system = TimeDependentQuantumSystem((a, t) -> H_drift(t) + a[1] * H_drive(t), n_drives)
-
-    @test system isa TimeDependentQuantumSystem
-    @test system.n_drives == n_drives
-    @test get_drift(system)(0.0) == H_drift(0.0)
-    @test get_drives(system)[1](0.0) == H_drive(0.0)
-
-    t = 2.0
-    @test get_drift(system)(t) == H_drift(t)
-    @test get_drives(system)[1](t) == H_drive(t)
-
-    # test three drives
-    n_drives = 3
-    Xₜ = t -> PAULIS.X * cos(t)
-    Yₜ = t -> PAULIS.Y * sin(t)
-    system = TimeDependentQuantumSystem((a, t) -> a[1] * Xₜ(t) + a[2] * Yₜ(t) + a[3] * PAULIS.Z, n_drives)
-
-    @test system isa TimeDependentQuantumSystem
-    @test get_drift(system)(0.0) == zeros(2, 2)
-    @test get_drift(system)(1.0) == zeros(2, 2)
-
-    t = 0.0
-    @test [H(t) for H in get_drives(system)] == [Xₜ(t), Yₜ(t), PAULIS.Z]
-
-    t = 2.0
-    @test [H(t) for H in get_drives(system)] == [Xₜ(t), Yₜ(t), PAULIS.Z]
-end
\ No newline at end of file
diff --git a/src/rollouts.jl b/src/rollouts.jl
index 3d18b43..5410384 100644
--- a/src/rollouts.jl
+++ b/src/rollouts.jl
@@ -187,10 +187,12 @@ function rollout(
 
     Ψ̃[:, 1] .= ψ̃_init
 
+    ts = cumsum([0.0; Δt[1:end-1]])
+
     p = Progress(T-1; enabled=show_progress)
     for t = 2:T
         aₜ₋₁ = controls[:, t - 1]
-        Gₜ = system.G(aₜ₋₁)
+        Gₜ = system.G(aₜ₋₁, ts[t - 1])
         if exp_vector_product
             Ψ̃[:, t] .= integrator(Δt[t - 1], Gₜ, Ψ̃[:, t - 1])
         else
@@ -488,10 +490,12 @@ function unitary_rollout(
 
     Ũ⃗[:, 1] .= Ũ⃗_init
 
+    ts = cumsum([0.0; Δt[1:end-1]])
+
     p = Progress(T-1; enabled=show_progress)
     for t = 2:T
         aₜ₋₁ = controls[:, t - 1]
-        Gₜ = system.G(aₜ₋₁)
+        Gₜ = system.G(aₜ₋₁, ts[t - 1])
         Ũₜ₋₁ = iso_vec_to_iso_operator(Ũ⃗[:, t - 1])
         if exp_vector_product
             Ũₜ = integrator(Δt[t - 1], Gₜ, Ũₜ₋₁)
@@ -916,7 +920,7 @@ end
     include("../test/test_utils.jl")
 
     traj = named_trajectory_type_1()
-    sys = QuantumSystem([PAULIS.X, PAULIS.Y])
+    sys = QuantumSystem([PAULIS.X, PAULIS.Y], 1.0, [(-1.0, 1.0), (-1.0, 1.0)])
     U_goal = GATES.H
     embedded_U_goal = EmbeddedOperator(U_goal, sys)
 
@@ -974,7 +978,7 @@ end
     using ForwardDiff
     using ExponentialAction
 
-    sys = QuantumSystem([PAULIS.X, PAULIS.Y])
+    sys = QuantumSystem([PAULIS.X, PAULIS.Y], 10.2, [(-1.0, 1.0), (-1.0, 1.0)])
     T = 51
     Δt = 0.2
     ts = fill(Δt, T)
@@ -991,7 +995,7 @@ end
     result2 = ForwardDiff.jacobian(
         as -> unitary_rollout(as, ts, sys, integrator=expv)[:, end], as
     )
-    iso_vec_dim = length(operator_to_iso_vec(sys.H(zeros(sys.n_drives))))
+    iso_vec_dim = length(operator_to_iso_vec(sys.H(zeros(sys.n_drives), 0.0)))
     @test size(result2) == (iso_vec_dim, T * sys.n_drives)
 
     # Time derivatives
@@ -1005,13 +1009,13 @@ end
     result2 = ForwardDiff.jacobian(
         ts -> unitary_rollout(as, ts, sys, integrator=expv)[:, end], ts
     )
-    iso_vec_dim = length(operator_to_iso_vec(sys.H(zeros(sys.n_drives))))
+    iso_vec_dim = length(operator_to_iso_vec(sys.H(zeros(sys.n_drives), 0.0)))
     @test size(result2) == (iso_vec_dim, T)
 end
 
 @testitem "Test variational rollouts" begin
     include("../test/test_utils.jl")
-    sys = QuantumSystem([PAULIS.X, PAULIS.Y])
+    sys = QuantumSystem([PAULIS.X, PAULIS.Y], 3.92, [(-1.0, 1.0), (-1.0, 1.0)])
     varsys1 = VariationalQuantumSystem([PAULIS.X, PAULIS.Y], [PAULIS.X])
     varsys2 = VariationalQuantumSystem([PAULIS.X, PAULIS.Y], [PAULIS.X, PAULIS.Y])
     U_goal = GATES.H
diff --git a/test_refactored_open_system.jl b/test_refactored_open_system.jl
new file mode 100644
index 0000000..e69de29