diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 27e73be..7e45265 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -19,15 +19,15 @@ jobs:
       matrix:
         version:
           - '1.6'
-          - '1.8'
+          - '1.10'
           - 'nightly'
         os:
           - ubuntu-latest
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
index 4230548..5336365 100644
--- a/.github/workflows/documentation.yml
+++ b/.github/workflows/documentation.yml
@@ -13,10 +13,10 @@ jobs:
       contents: write
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
-          version: '1.6'
+          version: '1.10'
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy
diff --git a/Project.toml b/Project.toml
index a9b9b83..d1917f4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,8 @@ version = "0.1.0"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+IntervalMatrices = "5c1f47dc-42dd-5697-8aaa-4d102d140ba9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RealTimeScheduling = "f1524149-62f5-490d-8249-9bfac76a7111"
@@ -15,6 +17,7 @@ RealTimeScheduling = "f1524149-62f5-490d-8249-9bfac76a7111"
 ControlSystemsBase = "1"
 DataStructures = "0.18"
 Distances = "0.10"
+IntervalArithmetic = "0.20"
 RealTimeScheduling = "0.3"
 julia = "^1.6"
 
diff --git a/README.md b/README.md
index e1dad2a..1c3ce11 100644
--- a/README.md
+++ b/README.md
@@ -54,8 +54,8 @@ notebook, `notebooks/control_timing_safety.jl`.
 ## Citing
 
 If you use this code in your research, we ask that you consider citing the
-relevant papers.  For the `bounded_runs_iter` function for overapproximating
-the maximum deviation, please cite:
+relevant papers.  For the `uncertain_linear_system` and `bounded_runs_iter`
+functions for overapproximating the maximum deviation, please cite:
 
 > Clara Hobbs, Bineet Ghosh, Shengjie Xu, Parasara Sridhar Duggirala, and
 > Samarjit Chakraborty.
diff --git a/docs/make.jl b/docs/make.jl
index d609191..5fa7f8c 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -7,6 +7,7 @@ makedocs(
     sitename = "ControlTimingSafety",
     format = Documenter.HTML(),
     modules = [ControlTimingSafety],
+    checkdocs = :exports,
     pages = [
         "index.md",
         "automata.md",
diff --git a/docs/src/safety.md b/docs/src/safety.md
index e54c649..7e4231e 100644
--- a/docs/src/safety.md
+++ b/docs/src/safety.md
@@ -1,13 +1,14 @@
 # Checking Safety
 
-The main safety-checking function of ControlTimingSafety is
-[`bounded_runs_iter`](@ref), which implements the Bounded Runs Iteration
-algorithm from "Safety Analysis of Embedded Controllers under Implementation
-Platform Timing Uncertainties."  This implementation permits an initial point,
-or an axis-aligned box initial set.  The helper function
+ControlTimingSafety implements two methods from "Safety Analysis of Embedded
+Controllers under Implementation Platform Timing Uncertainties": the Uncertain
+Linear Systems method in [`uncertain_linear_system`](@ref), and the Bounded
+Runs Iteration method in [`bounded_runs_iter`](@ref).  The implementation
+permits an axis-aligned box initial set.  The helper function
 [`bounded_runs`](@ref), which computes a single iteration, is also exported.
 
 ```@docs
+uncertain_linear_system
 bounded_runs_iter
 bounded_runs
 ```
diff --git a/notebooks/control_timing_safety.jl b/notebooks/control_timing_safety.jl
index 2449670..579100e 100644
--- a/notebooks/control_timing_safety.jl
+++ b/notebooks/control_timing_safety.jl
@@ -1,5 +1,5 @@
 ### A Pluto.jl notebook ###
-# v0.19.14
+# v0.19.22
 
 #> [frontmatter]
 #> title = "ControlTimingSafety.jl Demo"
@@ -29,12 +29,12 @@ begin
     Pkg.instantiate()
 
 	using Revise
-	using ControlTimingSafety
+	using ControlTimingSafety, IntervalArithmetic
 	using RealTimeScheduling
     using Plots, PlutoUI, LaTeXStrings
 	using Random, Distributions
 	using Distances
-	using ControlSystems, LinearAlgebra
+	using ControlSystemsBase, LinearAlgebra
 	using Printf
 end
 
@@ -150,7 +150,7 @@ let
 	for i = size(hit_prob, 1):-1:1
 	    seq = (rand(sim_time) .>= hit_prob[i]) .+ 1
 		z = evol(automaton, z_0, seq)
-	    plot!(z[:,1], z[:,2], label=L"P(\mathrm{hit}) = %$(hit_prob[i])", linecolor=line_colors[i])
+	    plot!(z[1,:], z[2,:], label=L"P(\mathrm{hit}) = %$(hit_prob[i])", linecolor=line_colors[i])
 	end
 	plot!()
 end
@@ -208,7 +208,7 @@ trj_2 = let
 	# Compute random trajectories
 	sp = SamplerUniformMissRow(MissRow(maxmiss_2), H_2)
 	seq = ones(Int64, H_2)
-	trj = Array{Float64}(undef, n_samples_2, H_2+1, size(T.Φ[1], 1))
+	trj = Array{Float64}(undef, n_samples_2, size(T.Φ[1], 1), H_2+1)
 	max_possible = (H_2 / (maxmiss_2 + 1)) * maxmiss_2 + H_2 % (maxmiss_2 + 1)
 	for i = axes(trj, 1)
 		accepted = false
@@ -232,7 +232,7 @@ We next compute the distance between each point of each random trajectory and th
 dist_2 = let	
 	dist = Array{Float64}(undef, n_samples_2, H_2+1)
 	for i in axes(trj_2, 1)
-		dist[i,:] = colwise(Euclidean(), trj_2[i,:,1:2]', nominal_2[:,1:2]')
+		dist[i,:] = colwise(Euclidean(), trj_2[i,1:2,:], nominal_2[1:2,:])
 	end
 	dist
 end
@@ -256,8 +256,8 @@ let
 	pipe_radius = sort(dev_2, dims=1)[end-1] - 0.01
 	circ_x = cos.(θ) * pipe_radius
 	circ_y = sin.(θ) * pipe_radius
-	for i in axes(nominal_2, 1)
-		plot!(circ_x .+ nominal_2[i,1], circ_y .+ nominal_2[i,2], repeat([i,], size(θ,1)), label=(i == 1) ? "Safety pipe" : "", seriestype=[:shape,], c=:lightblue, linecolor=:lightblue)
+	for i in axes(nominal_2, 2)
+		plot!(circ_x .+ nominal_2[1,i], circ_y .+ nominal_2[2,i], repeat([i,], size(θ,1)), label=(i == 1) ? "Safety pipe" : "", seriestype=[:shape,], c=:lightblue, linecolor=:lightblue)
 	end
 
 	# Plot random trajectories
@@ -265,7 +265,7 @@ let
 	# Plot the good trajectories first
 	for i = axes(trj_2, 1)
 		if dev_2[i] < pipe_radius
-			plot!(trj_2[i,:,1], trj_2[i,:,2], 1:H_2+1, label=(!label_good) ? "Random" : "", color=:green, opacity=10/n_samples_2)
+			plot!(trj_2[i,1,:], trj_2[i,2,:], 1:H_2+1, label=(!label_good) ? "Random" : "", color=:green, opacity=10/n_samples_2)
 			label_good = true
 		end
 	end
@@ -277,14 +277,14 @@ let
 			t_viol = argmax(dist_2[i,:], dims=1)
 			#label_bad || plot!(cos.(θ) * pipe_radius .+ nominal_2[t_viol,1], sin.(θ) * pipe_radius .+ nominal_2[t_viol,2], repeat([t_viol,], size(θ,1)), label="Radius violated", style=:dot, linecolor=:gray)
 			crosses = [(d >= pipe_radius) ? 1 : 0 for d in dist_2[i,:]]
-			plot!(trj_2[i,:,1], trj_2[i,:,2], 1:H_2+1, label=(!label_bad) ? "Violation" : "", color=:red, markershape=:x, markeralpha=crosses)
+			plot!(trj_2[i,1,:], trj_2[i,2,:], 1:H_2+1, label=(!label_bad) ? "Violation" : "", color=:red, markershape=:x, markeralpha=crosses)
 			nom_crosses .+= crosses
 			label_bad = true
 		end
 	end
 
 	# Finally, plot the nominal trajectory on top of everything else
-	plot!(nominal_2[:,1], nominal_2[:,2], 1:H_2+1, label="Nominal", color=:black, linewidth=3, marker=:x, markeralpha=nom_crosses)
+	plot!(nominal_2[1,:], nominal_2[2,:], 1:H_2+1, label="Nominal", color=:black, linewidth=3, marker=:x, markeralpha=nom_crosses)
 end
   ╠═╡ =#
 
@@ -314,7 +314,7 @@ This section illustrates the Bounded Runs algorithm, used to explore all possibl
 """
 
 # ╔═╡ 0843dd0e-0091-402f-9f28-fb112232f140
-bounds = [2; 2;; 2.5; 2.5]
+bounds = IntervalBox(2..2.5, 2..2.5)
 
 # ╔═╡ ef358d22-a530-4881-a282-b850d3655427
 begin
@@ -338,7 +338,7 @@ let
 	plot(title="Reachable set for $(n_bounded_runs) time steps", xlabel=L"x_1", ylabel=L"x_2", legend=:bottomright)
 	merged = merge_bounds(points_bounded_runs)
 	for k = 0:n_bounded_runs
-		corners = corners_from_bounds(merged[begin+k,:,:], cycle=true, dims=1:2)
+		corners = corners_from_bounds(merged[begin+k], cycle=true, dims=1:2)
 		plot!(corners[1,:], corners[2,:], label=L"x[%$(k)]")
 	end
 	
@@ -347,7 +347,7 @@ let
 		corners = augment(hsn, corners_from_bounds(bounds, dims=[1,2]))
 		for c in eachcol(corners)
 	    	x = evol(hsn, c, ones(Int64, n_bounded_runs))
-			plot!(x[:,1], x[:,2], label="Sample", linecolor=:blue, marker=:circle)
+			plot!(x[1,:], x[2,:], label="Sample", linecolor=:blue, marker=:circle)
 		end
 	end
 	plot!()
@@ -390,7 +390,7 @@ end
 let
 	plot(title="Reachable set for $(time_4) time steps, $(n_4) per tree", xlabel=L"x_1", ylabel=L"x_2", legend=:bottomright)
 	for t in 1:time_4+1
-		corners = corners_from_bounds(all_bounds[t,:,:], cycle=true, dims=1:2)
+		corners = corners_from_bounds(all_bounds[t], cycle=true, dims=1:2)
 		plot!(corners[1,:], corners[2,:], label=L"x[%$(t-1)]")
 	end
 	
@@ -399,7 +399,7 @@ let
 		corners = augment(hsn, corners_from_bounds(bounds, dims=[1,2]))
 		for c in eachcol(corners)
 	    	x = evol(hsn, c, ones(Int64, n_4*t_4))
-			plot!(x[:,1], x[:,2], label="Sample", linecolor=:blue, marker=:circle)
+			plot!(x[1,:], x[2,:], label="Sample", linecolor=:blue, marker=:circle)
 		end
 	end
 	plot!(legend=false)
@@ -413,7 +413,7 @@ Using these computed reachable sets, we finally compute the maximum deviation th
 # ╔═╡ 7afec811-9ccc-4e93-9cc8-a3c6b6ee5c49
 let
 	automaton = strat_map[sim_strategy_4](sysd, K, MissRow(max_miss_4))
-	d = deviation(automaton, augment(automaton, bounds), all_bounds, dims=1:2)[1:time_4+1]
+	d = deviation(automaton, augment(automaton, bounds), all_bounds)[1:time_4+1]
 	i = argmax(d)
 	v = maximum(d)
 	
diff --git a/src/ControlTimingSafety.jl b/src/ControlTimingSafety.jl
index 6f0bb99..058ade1 100644
--- a/src/ControlTimingSafety.jl
+++ b/src/ControlTimingSafety.jl
@@ -4,6 +4,7 @@ using LinearAlgebra
 using Random
 using ControlSystemsBase
 using Distances
+using IntervalArithmetic
 
 export Automaton, nlocations, nactions
 export hold_skip_next, zero_skip_next, hold_kill, zero_kill
@@ -11,6 +12,7 @@ export strat_map, strat_names
 export evol_final, evol, evol_final!, evol!, augment
 include("automata.jl")
 
+export uncertain_linear_system
 export corners_from_bounds, merge_bounds, merge_bounds!
 export bounded_runs, bounded_runs_iter, deviation
 include("safety.jl")
diff --git a/src/automata.jl b/src/automata.jl
index 7ece6a3..72be85a 100644
--- a/src/automata.jl
+++ b/src/automata.jl
@@ -397,19 +397,47 @@ and `l` is the final location in the automaton.
 """
 function evol_final(a::Automaton, z_0::AbstractVector{Float64}, input::AbstractVector{Int64})
     @boundscheck length(z_0) == a.nz || throw(DimensionMismatch("z_0 must have length a.nz"))
-    z = zeros(length(input) + 1, a.nz)
-    z[1,:] = z_0
+    z = zeros(a.nz, length(input) + 1)
+    z[:,1] = z_0
+    evol_final!(a, z, input)
+end
+
+function evol_final(a::Automaton, z_0::IntervalBox, input::AbstractVector{Int64})
+    @boundscheck length(z_0) == a.nz || throw(DimensionMismatch("z_0 must have length a.nz"))
+    z = Vector{IntervalBox{length(z_0), eltype(inf.(z_0))}}(undef, length(input) + 1)
+    z[1] = z_0
     evol_final!(a, z, input)
 end
 
 """
     evol_final!(a, z, input)
 
-Same as [`evol_final`](@ref), but writes to the input matrix `z`, whose first row is `z_0`.
+Same as [`evol_final`](@ref), but writes to the input matrix `z`, whose first column is `z_0`.
+`z` may also be a vector of `IntervalBox`es, whose first element is `z_0`.
 """
 function evol_final!(a::Automaton, z::AbstractMatrix{Float64}, input::AbstractVector{Int64})
-    @boundscheck size(z,1) == length(input)+1 || throw(DimensionMismatch("z must have size (length(input)+1, a.nz)"))
-    @boundscheck size(z,2) == a.nz || throw(DimensionMismatch("z must have size (length(input)+1, a.nz)"))
+    @boundscheck size(z,1) == a.nz || throw(DimensionMismatch("z must have size (a.nz, length(input)+1)"))
+    @boundscheck size(z,2) == length(input)+1 || throw(DimensionMismatch("z must have size (a.nz, length(input)+1)"))
+    l = a.l_int
+    # For each time step and input action
+    for (t, in) in enumerate(input)
+        # Get the dynamics matrix
+        μ = a.μ[l, in]
+        # If we hit a missing transition, return the states that we reached,
+        # and a missing final location to signal the problem to the caller.
+        if ismissing(μ)
+            return z[:,1:t], missing
+        end
+        # Apply the dynamics
+        mul!(view(z,:,t+1), a.Φ[μ], view(z,:,t))
+        # Transition to the new location
+        l = a.T[l, in]
+    end
+    z, l
+end
+
+function evol_final!(a::Automaton, z::AbstractVector{IntervalBox{S, T}}, input::AbstractVector{Int64}) where {S, T}
+    @boundscheck length(z) == length(input)+1 || throw(DimensionMismatch("z must have length length(input)+1"))
     l = a.l_int
     # For each time step and input action
     for (t, in) in enumerate(input)
@@ -418,10 +446,10 @@ function evol_final!(a::Automaton, z::AbstractMatrix{Float64}, input::AbstractVe
         # If we hit a missing transition, return the states that we reached,
         # and a missing final location to signal the problem to the caller.
         if ismissing(μ)
-            return z[1:t,:], missing
+            return z[1:t], missing
         end
         # Apply the dynamics
-        mul!(view(z,t+1,:), a.Φ[μ], view(z,t,:))
+        z[t+1] = a.Φ[μ] * z[t]
         # Transition to the new location
         l = a.T[l, in]
     end
@@ -439,14 +467,14 @@ Returns `z`, a matrix of states over time.
 See also [`evol_final`](@ref), which additionally returns the final location in the
 automaton.
 """
-evol(a::Automaton, z_0::AbstractVector{Float64}, input::AbstractVector{Int64}) = evol_final(a, z_0, input)[1]
+evol(a::Automaton, z_0::Union{AbstractVector{Float64}, IntervalBox{S, T}}, input::AbstractVector{Int64}) where {S, T} = evol_final(a, z_0, input)[1]
 
 """
     evol!(a, z, input)
 
-Same as [`evol`](@ref), but writes to the input matrix `z`, whose first row is `z_0`.
+Same as [`evol`](@ref), but writes to the input matrix `z`, whose first column is `z_0`.
 """
-evol!(a::Automaton, z::AbstractMatrix{Float64}, input::AbstractVector{Int64}) = evol_final!(a, z, input)[1]
+evol!(a::Automaton, z::Union{AbstractMatrix{Float64}, AbstractVector{IntervalBox{S, T}}}, input::AbstractVector{Int64}) where {S, T} = evol_final!(a, z, input)[1]
 
 """
     augment(a::Automaton, x)
@@ -456,3 +484,4 @@ augmented state space of `a`.
 """
 augment(a::Automaton, x::AbstractMatrix) = [x; zeros(a.nz - size(x, 1), size(x, 2))]
 augment(a::Automaton, x::AbstractVector) = [x; zeros(a.nz - size(x, 1))]
+augment(a::Automaton, x::IntervalBox) = IntervalBox([x...; zeros(a.nz - length(x))])
diff --git a/src/probablesafety.jl b/src/probablesafety.jl
index 00d8c97..4a85d9e 100644
--- a/src/probablesafety.jl
+++ b/src/probablesafety.jl
@@ -9,19 +9,18 @@ as in [`deviation`](@ref).  If `estimate` is specified, stop early if this devia
 estimate is exceeded, returning the exceeding deviation.
 """
 Base.@propagate_inbounds function maximum_deviation_random(a::Automaton,
-        sampler::RealTimeScheduling.SamplerWeaklyHard, samples::Int64,
-        z_0::AbstractVecOrMat{Float64};
+        sampler::RealTimeScheduling.SamplerWeaklyHard, samples::Int64, z_0::IntervalBox;
         estimate::Union{Float64, Nothing}=nothing, metric::PreMetric=Euclidean(),
         nominal::AbstractVector{Int64}=ones(Int64,sampler.H),
         nominal_trajectory::Union{AbstractArray{Float64}, Nothing}=nothing)
     @boundscheck length(nominal) == sampler.H || throw(DimensionMismatch("nominal must have length sampler.H"))
-    @boundscheck nominal_trajectory === nothing || size(nominal_trajectory) == (size(z_0,1), sampler.H+1) || throw(DimensionMismatch("nominal_trajectory must have size (size(z_0,1), sampler.H+1)"))
+    @boundscheck nominal_trajectory === nothing || size(nominal_trajectory) == (a.nz, sampler.H+1) || throw(DimensionMismatch("nominal_trajectory must have size (size(z_0,1), sampler.H+1)"))
     corners = unique(corners_from_bounds(z_0), dims=2)
     # Compute the nominal trajectory if not provided
     if nominal_trajectory === nothing
-        nominal_trajectory = Array{Float64}(undef, size(z_0,1), sampler.H+1, size(corners,2))
+        nominal_trajectory = Array{Float64}(undef, a.nz, sampler.H+1, size(corners,2))
         for i in axes(corners, 2)
-            nominal_trajectory[:,:,i] = evol(a, corners[:,i], nominal)'
+            nominal_trajectory[:,:,i] = evol(a, corners[:,i], nominal)
         end
     end
     Cnom = Array{Float64}(undef, size(a.C,1), sampler.H+1, size(corners,2))
@@ -33,13 +32,13 @@ Base.@propagate_inbounds function maximum_deviation_random(a::Automaton,
     s = falses(sampler.H)
     beh = zeros(Int64, sampler.H)
     if z_0 isa AbstractVector
-        z = zeros(sampler.H + 1, a.nz, 1)
+        z = zeros(a.nz, sampler.H + 1, 1)
         Cz = zeros(size(a.C, 1), sampler.H + 1, 1)
     else
-        z = zeros(sampler.H + 1, a.nz, size(corners, 2))
+        z = zeros(a.nz, sampler.H + 1, size(corners, 2))
         Cz = zeros(size(a.C, 1), sampler.H + 1, size(z, 3))
     end
-    z[1,:,:] = corners
+    z[:,1,:] = corners
     dist = zeros(Float64, length(nominal)+1)
 
     H = zeros(sampler.H+1)
@@ -54,7 +53,7 @@ Base.@propagate_inbounds function maximum_deviation_random(a::Automaton,
             # Compute its evolution
             evol!(a, view(z,:,:,i), beh)
             # Compute the output
-            mul!(view(Cz,:,:,i), a.C, view(z,:,:,i)')
+            mul!(view(Cz,:,:,i), a.C, view(z,:,:,i))
         end
         # Compute the distance between the output and nominal trajectory
         for t in eachindex(H)
@@ -71,6 +70,19 @@ Base.@propagate_inbounds function maximum_deviation_random(a::Automaton,
     maximum(H)
 end
 
+Base.@propagate_inbounds function maximum_deviation_random(a::Automaton,
+        sampler::RealTimeScheduling.SamplerWeaklyHard, samples::Int64, z_0::Vector;
+        estimate::Union{Float64, Nothing}=nothing, metric::PreMetric=Euclidean(),
+        nominal::AbstractVector{Int64}=ones(Int64,sampler.H),
+        nominal_trajectory::Union{AbstractArray{Float64}, Nothing}=nothing)
+    maximum_deviation_random(a, sampler, samples, IntervalBox(z_0);
+            estimate=estimate,
+            metric=metric,
+            nominal=nominal,
+            nominal_trajectory=nominal_trajectory)
+end
+
+
 """
     estimate_deviation(a, sampler, z_0, c, B; nominal, nominal_trajectory,
                        estimate, estimate_samples, ϵ)
diff --git a/src/safety.jl b/src/safety.jl
index e3e624d..e6f0e8b 100644
--- a/src/safety.jl
+++ b/src/safety.jl
@@ -1,19 +1,46 @@
+using IntervalArithmetic, IntervalMatrices
+
+
 """
-    corners_from_bounds(bounds::AbstractMatrix; cycle=false, dims=axes(bounds, 1))
+    uncertain_linear_system(a::Automaton, z_0::IntervalBox, t)
 
-Returns the corners of the n-dimensional interval represented by `bounds`.  If `cycle` is
-`true`, the first corner is repeated at the end, and the corners are given in Gray code
-order.  Only the dimensions from `dims` are considered.
+Compute reachable sets for `t` time steps for the given [`Automaton`](@ref) `a`, starting
+from the initial set given by `z_0`.  `z_0` must have length `a.nz`.
 """
-function corners_from_bounds(bounds::AbstractMatrix; cycle::Bool=false, dims=axes(bounds, 1))
-    @boundscheck dims ⊆ axes(bounds, 1) || throw(ArgumentError("All entries of dims must be valid indices to the first dimension of bounds"))
-    if size(bounds, 2) == 1
-        return bounds
+function uncertain_linear_system(a::Automaton, z_0::IntervalBox, t)
+    # Compute the uncertain dynamics matrix Λ
+    Φ = cat(a.Φ..., dims=3)
+    mins = minimum(Φ, dims=3)
+    maxs = maximum(Φ, dims=3)
+    shape = (size(mins, 1), size(mins, 2))
+    Λ = IntervalMatrix([m..M for (m, M) in zip(reshape(mins, shape), reshape(maxs, shape))])
+
+    # Compute the reachable sets
+    all_bounds = Array{IntervalBox}(undef, t+1)
+    all_bounds[1] = z_0
+    for i in eachindex(all_bounds)
+        i > 1 || continue
+        all_bounds[i] = Λ * all_bounds[i-1]
     end
+
+    all_bounds
+end
+
+
+"""
+    corners_from_bounds(bounds::IntervalBox; cycle=false, dims=1:length(bounds))
+
+Returns the corners of the n-dimensional interval `bounds`.  If `cycle` is `true`, the first
+corner is repeated at the end, and the corners are given in Gray code order.  Only the
+dimensions from `dims` are considered.
+"""
+function corners_from_bounds(bounds::IntervalBox; cycle::Bool=false, dims=Base.oneto(length(bounds)))
+    @boundscheck dims ⊆ 1:length(bounds) || throw(ArgumentError("All entries of dims must be valid indices to the first dimension of bounds"))
     ldims = length(dims)
 
-    corners = cat(reshape([[c...] for c in Base.product(eachrow(bounds[dims,:])...)],
-              2^ldims)..., dims=2)
+    corners = cat(reshape([[c...] for c in Base.product(
+                           eachrow([inf.(bounds[dims]) sup.(bounds[dims])])...)],
+                  2^ldims)..., dims=2)
     if cycle
         gray(x) = x ⊻ (x >> 1)
         [corners[:,gray.(0:2^ldims-1) .+ 1] corners[:,1]]
@@ -25,12 +52,12 @@ end
 """
     corners_from_bounds(bounds::AbstractVector; cycle=nothing, dims=nothing)
 
-When applied to a vector, cast it to a one-column matrix.  `cycle` and `dims` are ignored.
+When applied to a vector, returns a 3-dimensional `Array` whose third dimension corresponds
+to the index of the input `Vector`.
 """
-corners_from_bounds(bounds::AbstractVector; cycle=nothing, dims=nothing) = reshape(bounds, length(bounds), 1)
-
-_safefloatmin(x) = (isnan(x) || isinf(x)) ? +Inf : x
-_safefloatmax(x) = (isnan(x) || isinf(x)) ? -Inf : x
+function corners_from_bounds(bounds::AbstractVector; cycle::Bool=false, dims=Base.oneto(length(bounds[begin])))
+    cat(corners_from_bounds.(bounds, cycle=cycle, dims=dims)..., dims=3)
+end
 
 """
     merge_bounds(b)
@@ -40,7 +67,7 @@ Merges an array of bounding boxes `b` into one.
 See also [`merge_bounds!`](@ref).
 """
 function merge_bounds(b)
-    r = similar(b, axes(b)[2:4])
+    r = similar(b, axes(b, 1))
     merge_bounds!(r, b)
 end
 
@@ -53,24 +80,23 @@ Merges an array of bounding boxes `b` into one, storing the result in `r`.
 See also [`merge_bounds`](@ref).
 """
 function merge_bounds!(r, b)
-    minimum!(_safefloatmin, reshape(view(r, :, :, 1), 1, size(r, 1), size(r, 2)), view(b, :, :, :, 1))
-    maximum!(_safefloatmax, reshape(view(r, :, :, 2), 1, size(r, 1), size(r, 2)), view(b, :, :, :, 2))
+    for t in eachindex(r)
+        r[t] = hull(b[t,:]...)
+    end
     r
 end
 
 mutable struct _StackFrame
-    z::Matrix{Float64}
+    z::IntervalBox
     loc::Int64
     act::Int64
 end
+
 """
-    bounded_runs(a::Automaton, z_0, n)
+    bounded_runs(a::Automaton, z_0::IntervalBox, n)
 
 Compute reachable sets for `n` time steps for the given [`Automaton`](@ref) `a`, starting
-from the initial set given by `z_0`.
-
-`z_0` must be an `a.nz`-element vector, or an `a.nz`×`2` matrix whose first and second
-columns are the minimum and maximum along each dimension, respectively.
+from the initial set given by `z_0`.  `z_0` must have length `a.nz`.
 
 Returns `(bounds, locs)`, where `bounds` is an `nactions(a)^n`×`n+1`×`a.nz`×`2`
 `Array{Float64}` giving the bounding box for each (run, time step), and `locs` is an
@@ -82,23 +108,19 @@ compute reachable sets for longer time horizons.  Typically one will call
 [`deviation`](@ref) on the results of this function to determine deviation from a nominal
 trajectory.
 """
-function bounded_runs(a::Automaton, z_0::AbstractVecOrMat, n::Integer)
-    corners = corners_from_bounds(z_0)
-
+function bounded_runs(a::Automaton, z_0::IntervalBox, n::Integer)
     # Stack
-    # z gets one extra entry in the third dimension for cheap concatenation in leaf nodes
     st = Array{_StackFrame}(undef, n+1)
     for i in eachindex(st)
-        st[i] = _StackFrame(Array{Float64}(undef, size(corners,1), size(corners,2)+1), a.l_int, 1)
+        st[i] = _StackFrame(IntervalBox(∅, a.nz), a.l_int, 1)
     end
 
     # Bounding boxes for each time step, final location
-    ret = Array{Float64}(undef, a.nz, 2, n+1, nlocations(a))
-    ret[:,1,:,:] .= Inf
-    ret[:,2,:,:] .= -Inf
+    ret = Array{IntervalBox}(undef, n+1, nlocations(a))
+    ret[:,:] .= repeat([IntervalBox(∅, a.nz)], n+1, nlocations(a))
 
     # Create the stack frame for time 0
-    st[1].z[:,begin:end-1] = corners
+    st[1].z = z_0
     # Initialize the stack pointer
     sp = 1
     # While we haven't popped all the way out
@@ -107,10 +129,7 @@ function bounded_runs(a::Automaton, z_0::AbstractVecOrMat, n::Integer)
         if sp == n+1
             # Calculate min and max for this final location at each time step
             for (i, sf) in enumerate(st)
-                sf.z[:,end] = view(ret,:,1,i,st[sp].loc)
-                minimum!(view(ret,:,1:1,i,st[sp].loc), sf.z, init=false)
-                sf.z[:,end] = view(ret,:,2,i,st[sp].loc)
-                maximum!(view(ret,:,2:2,i,st[sp].loc), sf.z, init=false)
+                ret[i,st[sp].loc] = hull(ret[i, st[sp].loc], sf.z)
             end
             sp -= 1
         # If we're out of actions from this step
@@ -122,20 +141,18 @@ function bounded_runs(a::Automaton, z_0::AbstractVecOrMat, n::Integer)
             st[sp].act += 1
         # If the transition is present
         else
-            mul!(st[sp+1].z, a.Φ[a.μ[st[sp].loc, st[sp].act]], st[sp].z)
+            st[sp+1].z = a.Φ[a.μ[st[sp].loc, st[sp].act]] * st[sp].z
             st[sp+1].loc = a.T[st[sp].loc, st[sp].act]
             st[sp+1].act = 1
             st[sp].act += 1
             sp += 1
         end
     end
-    # TODO: the order of dimensions is largely an implementation detail, and since it was
-    # non-optimal before, we should ultimately remove the need for this permutedims call
-    permutedims(ret, [4, 3, 1, 2])
+    ret
 end
 
 """
-    bounded_runs_iter(a, z_0, n, t)
+    bounded_runs_iter(a::Automaton, z_0::IntervalBox, n, t)
 
 Iterate [`bounded_runs`](@ref)`(a, z_0, n)` for `t` iterations, returning the reachable
 set at each of the `n`×`t+1` time steps.
@@ -143,17 +160,13 @@ set at each of the `n`×`t+1` time steps.
 See also [`deviation`](@ref), which can be called with the result of this function to find
 the deviation from a nominal trajectory.
 """
-function bounded_runs_iter(a::Automaton, z_0::AbstractVecOrMat, n::Integer, t::Integer; safety_margin::Float64=Inf)
-    # Dimensions: time, augmented state, min/max
-    all_bounds = Array{Float64}(undef, n*(t+1)+1, a.nz, 2)
-    if isa(z_0, AbstractVector)
-        all_bounds[1,:,:] = [z_0 z_0]
-    else
-        all_bounds[1,:,:] = z_0
-    end
+function bounded_runs_iter(a::Automaton, z_0::IntervalBox, n::Integer, t::Integer; safety_margin::Float64=Inf)
+    # Dimensions: time
+    all_bounds = Array{IntervalBox}(undef, n*(t+1)+1)
+    all_bounds[1] = z_0
 
     bounds = bounded_runs(a, z_0, n)
-    merge_bounds!(view(all_bounds, 1:n+1, :, :), bounds)
+    merge_bounds!(view(all_bounds, 1:n+1), bounds)
 
     A = Array{Automaton}(undef, length(a.L))
     for i in a.L
@@ -162,36 +175,36 @@ function bounded_runs_iter(a::Automaton, z_0::AbstractVecOrMat, n::Integer, t::I
 
     if isfinite(safety_margin)
         nominal = ones(Int64, n*(t+1))
-        nom = Array{Float64}(undef, size(a.C, 1), 2^size(z_0,1), n*(t+1)+1)
+        nom = Array{Float64}(undef, size(a.C, 1), 2^size(a.C,2), n*(t+1)+1)
         corners = corners_from_bounds(z_0)
         for (i, c) in enumerate(eachcol(corners))
             e = evol(a, c, nominal)
-            nom[:,i,:] = a.C * e'
+            nom[:,i,:] = a.C * e
         end
-        d = deviation(a, z_0, all_bounds[1:n+1,:,:], nominal_trajectory=nom[:,:,1:n+1])
+        d = deviation(a, z_0, all_bounds[1:n+1], nominal_trajectory=nom[:,:,1:n+1])
         if maximum(d) > safety_margin
-            return all_bounds[1:n+1,:,:]
+            return all_bounds[1:n+1]
         end
     end
 
-    # Dimensions: initial location, final location, time, augmented state, min/max
-    new_bounds = Array{Float64}(undef, nlocations(a), nlocations(a), n+1, a.nz, 2)
+    # Dimensions: initial location, final location, time
+    new_bounds = Array{IntervalBox}(undef, n+1, nlocations(a), nlocations(a))
     for i in 1:t
         # Simulate each box from previous iteration
         Threads.@threads for i in a.L
-            new_bounds[i,:,:,:,:] = bounded_runs(A[i], bounds[i,end,:,:], n)
+            new_bounds[:,:,i] = bounded_runs(A[i], bounds[end,i], n)
         end
         # Merge resulting boxes from these simulations
         for i in a.L
-            merge_bounds!(view(bounds, i, :, :, :), view(new_bounds, :, i, :, :, :))
+            merge_bounds!(view(bounds, :, i), view(new_bounds, :, i, :))
         end
         # Save the bounds
-        merge_bounds!(view(all_bounds, n*i+1:n*(i+1)+1, :, :), bounds)
+        merge_bounds!(view(all_bounds, n*i+1:n*(i+1)+1), bounds)
 
         if isfinite(safety_margin)
-            d = deviation(a, z_0, all_bounds[n*i+2:n*(i+1)+1,:,:], nominal_trajectory=nom[:,:,n*i+2:n*(i+1)+1])
+            d = deviation(a, z_0, all_bounds[n*i+2:n*(i+1)+1], nominal_trajectory=nom[:,:,n*i+2:n*(i+1)+1])
             if maximum(d) > safety_margin
-                return all_bounds[1:n*(i+1)+1,:,:]
+                return all_bounds[1:n*(i+1)+1]
             end
         end
     end
@@ -216,27 +229,23 @@ trajectory from the `nominal` behavior.  This can improve efficiency when callin
 See also [`bounded_runs`](@ref) and [`bounded_runs_iter`](@ref), which can be used to
 compute `reachable`.
 """
-Base.@propagate_inbounds function deviation(a::Automaton, z_0::AbstractVecOrMat{Float64},
-                   reachable::AbstractArray{Float64,3};
+Base.@propagate_inbounds function deviation(a::Automaton, z_0::IntervalBox,
+                   reachable::AbstractVector{IntervalBox};
                    metric::PreMetric=Euclidean(),
                    nominal::AbstractVector{Int64}=ones(Int64,size(reachable,1)-1),
                    nominal_trajectory::Union{AbstractArray{Float64,3}, Nothing}=nothing)
-    @boundscheck length(nominal) == size(reachable, 1) - 1 || throw(DimensionMismatch("nominal must have length size(reachable, 1) - 1"))
+    @boundscheck length(nominal) == length(reachable) - 1 || throw(DimensionMismatch("nominal must have length length(reachable) - 1"))
 
     # Dimensions: state variables, points, time
-    reachable_corners = cat([corners_from_bounds(a.C * reachable[t,:,:]) for t in axes(reachable, 1)]..., dims=3)
+    reachable_corners = corners_from_bounds([a.C * reachable[t] for t in eachindex(reachable)])
 
     # Dimensions: state variables, points, time
     if nominal_trajectory === nothing
-        if z_0 isa AbstractVector{Float64}
-            nominal_trajectory = reshape(a.C * evol(a, z_0, nominal)', size(a.C, 1), 1, size(reachable, 1))
-        else
-            nominal_trajectory = Array{Float64}(undef, size(a.C, 1), 2^size(z_0,1), size(reachable, 1))
-            corners = corners_from_bounds(z_0)
-            for (i, c) in enumerate(eachcol(corners))
-                e = evol(a, c, nominal)
-                nominal_trajectory[:,i,:] = a.C * e'
-            end
+        nominal_trajectory = Array{Float64}(undef, size(a.C, 1), 2^length(z_0), length(reachable))
+        corners = corners_from_bounds(z_0)
+        for (i, c) in enumerate(eachcol(corners))
+            e = evol(a, c, nominal)
+            nominal_trajectory[:,i,:] = a.C * e
         end
     end
 
diff --git a/src/schedule_synthesis.jl b/src/schedule_synthesis.jl
index ce5522f..256200a 100644
--- a/src/schedule_synthesis.jl
+++ b/src/schedule_synthesis.jl
@@ -80,7 +80,7 @@ guarantees the deviation upper bound is at most `d_max`. The system is specified
 [`bounded_runs_iter`](@ref).
 """
 function synthesize_constraints(sysd::AbstractStateSpace{<:Discrete},
-    K::AbstractMatrix{Float64}, z_0::AbstractVecOrMat, d_max::Float64,
+    K::AbstractMatrix{Float64}, z_0::IntervalBox, d_max::Float64,
     maxwindow::Int64, n::Int64, t::Int64)
 
     safe_constraints = MeetAny[]
diff --git a/test/Project.toml b/test/Project.toml
index 949cea8..324ec53 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
+IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 RealTimeScheduling = "f1524149-62f5-490d-8249-9bfac76a7111"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/schedule_synthesis_tests.jl b/test/schedule_synthesis_tests.jl
index 783c481..7166bea 100644
--- a/test/schedule_synthesis_tests.jl
+++ b/test/schedule_synthesis_tests.jl
@@ -1,5 +1,6 @@
 using LinearAlgebra
 using RealTimeScheduling
+using IntervalArithmetic
 
 @testset "Schedule Synthesis" begin
      sys = let
@@ -35,7 +36,7 @@ using RealTimeScheduling
           lqr(sysd_delay, Q, R)
      end
 
-     @test synthesize_constraints(sysd, K, repeat([10.0], 3), 5.5, 6, 10, 10) == [
+     @test synthesize_constraints(sysd, K, IntervalBox(10..10, 10..10, 10..10), 5.5, 6, 10, 10) == [
           RealTimeScheduling.MeetAny(1, 2),
           RealTimeScheduling.MeetAny(2, 3),
           RealTimeScheduling.MeetAny(3, 4),