diff --git a/lib/DiffEqDevTools/src/DiffEqDevTools.jl b/lib/DiffEqDevTools/src/DiffEqDevTools.jl
index bc43670fda9..5dd0696ce46 100644
--- a/lib/DiffEqDevTools/src/DiffEqDevTools.jl
+++ b/lib/DiffEqDevTools/src/DiffEqDevTools.jl
@@ -32,6 +32,7 @@ const ALL_ERRORS = union([:final],
 include("benchmark.jl")
 include("convergence.jl")
 include("plotrecipes.jl")
+include("autoplot.jl")
 include("test_solution.jl")
 include("ode_tableaus.jl")
 include("tableau_info.jl")
@@ -45,6 +46,21 @@ export test_convergence, analyticless_test_convergence, appxtrue!, appxtrue
 
 export get_sample_errors
 
+#Tagging and Filtering
+export filter_by_tags, exclude_by_tags, get_tags, unique_tags, merge_wp_sets
+
+#Multi-error-mode
+export available_errors
+
+#Best-of-family helpers
+export best_by_tag, best_of_families, wp_area
+
+#AutoDiff comparison helpers
+export with_autodiff_variants
+
+#Autoplot
+export autoplot
+
 #Tab Functions
 export stability_region, residual_order_condition, check_tableau, imaginary_stability_interval
 
diff --git a/lib/DiffEqDevTools/src/autoplot.jl b/lib/DiffEqDevTools/src/autoplot.jl
new file mode 100644
index 00000000000..2129d52c678
--- /dev/null
+++ b/lib/DiffEqDevTools/src/autoplot.jl
@@ -0,0 +1,81 @@
+"""
+    autoplot(wp_set; families=nothing, reference_tags=nothing, best_n=2)
+
+Generate a collection of WorkPrecisionSets for standard comparison plots.
+
+Returns a `Dict{String, WorkPrecisionSet}` with:
+- `"family_<tag>"` for each family
+- `"best_of_families"` with top methods per family
+- `"all"` containing the original set
+"""
+function autoplot(wp_set::WorkPrecisionSet;
+        families = nothing,
+        reference_tags = nothing,
+        best_n::Int = 2)
+    results = Dict{String, WorkPrecisionSet}()
+
+    # Determine families
+    if families === nothing
+        families = _auto_detect_families(wp_set)
+    end
+
+    # Per-family plots
+    for fam in families
+        subset = filter_by_tags(wp_set, fam)
+        if reference_tags !== nothing
+            ref_tags = reference_tags isa Symbol ? (reference_tags,) : Tuple(reference_tags)
+            ref = filter_by_tags(wp_set, ref_tags...)
+            # Deduplicate
+            seen = Set(wp.name for wp in subset.wps)
+            extra_wps = [wp for wp in ref.wps if wp.name ∉ seen]
+            if !isempty(extra_wps)
+                subset = merge_wp_sets(subset,
+                    WorkPrecisionSet(extra_wps, length(extra_wps),
+                        wp_set.abstols, wp_set.reltols, wp_set.prob,
+                        wp_set.setups, [wp.name for wp in extra_wps],
+                        wp_set.error_estimate, wp_set.numruns))
+            end
+        end
+        length(subset) > 0 && (results["family_$(fam)"] = subset)
+    end
+
+    # Best-of-families
+    if !isempty(families)
+        best_wps = WorkPrecision[]
+        seen = Set{String}()
+        for fam in families
+            subset = filter_by_tags(wp_set, fam)
+            length(subset) == 0 && continue
+            ranked = sort(collect(1:length(subset.wps)),
+                by = i -> let t = filter(!isnan, subset.wps[i].times)
+                    isempty(t) ? Inf : mean(t)
+                end)
+            for idx in ranked[1:min(best_n, length(ranked))]
+                wp = subset.wps[idx]
+                if wp.name ∉ seen
+                    push!(seen, wp.name)
+                    push!(best_wps, wp)
+                end
+            end
+        end
+        if !isempty(best_wps)
+            names = [wp.name for wp in best_wps]
+            results["best_of_families"] = WorkPrecisionSet(
+                best_wps, length(best_wps), wp_set.abstols, wp_set.reltols,
+                wp_set.prob, wp_set.setups, names, wp_set.error_estimate, wp_set.numruns)
+        end
+    end
+
+    results["all"] = wp_set
+    return results
+end
+
+function _auto_detect_families(wp_set::WorkPrecisionSet)
+    all_tags = unique_tags(wp_set)
+    n_methods = length(wp_set.wps)
+    n_methods == 0 && return Symbol[]
+    # Filter out tags that appear on >80% of methods (likely category tags)
+    return filter(all_tags) do tag
+        count(wp -> tag in wp.tags, wp_set.wps) <= 0.8 * n_methods
+    end
+end
diff --git a/lib/DiffEqDevTools/src/benchmark.jl b/lib/DiffEqDevTools/src/benchmark.jl
index 0babc076815..69bb1673385 100644
--- a/lib/DiffEqDevTools/src/benchmark.jl
+++ b/lib/DiffEqDevTools/src/benchmark.jl
@@ -1,5 +1,17 @@
 using Statistics
 
+# Convert a vector of error Dicts to a StructArray with consistent NamedTuple key ordering.
+# Dict iteration order is not guaranteed, so NamedTuple.(dicts) can produce NamedTuples
+# with different type parameters, causing the element type to widen to bare `NamedTuple`.
+# StructArrays 0.7 requires parameterized NamedTuple types (NamedTuple{names}), so we
+# must ensure all NamedTuples share the same key ordering.
+function _dicts_to_structarray(dicts)
+    ks = Tuple(sort(collect(keys(first(dicts)))))
+    V = valtype(first(dicts))
+    NT = NamedTuple{ks, NTuple{length(ks), V}}
+    return StructArray([NT(Tuple(d[k] for k in ks)) for d in dicts])
+end
+
 # Default names of algorithms:
 # Workaround for `MethodOfSteps` algorithms, otherwise they are all called "MethodOfSteps"
 # Ideally this would be a trait (in SciMLBase?), so packages could implement it
@@ -167,6 +179,7 @@ mutable struct WorkPrecision
     name::Any
     error_estimate::Any
     N::Int
+    tags::Vector{Symbol}
 end
 
 mutable struct WorkPrecisionSet
@@ -179,11 +192,37 @@ mutable struct WorkPrecisionSet
     names::Any
     error_estimate::Any
     numruns::Any
+    active_error_estimates::Vector{Symbol}
+end
+
+# Backward-compatible 9-argument constructor
+function WorkPrecisionSet(wps, N, abstols, reltols, prob, setups, names, error_estimate, numruns)
+    WorkPrecisionSet(wps, N, abstols, reltols, prob, setups, names, error_estimate, numruns,
+        [error_estimate isa Symbol ? error_estimate : :final])
+end
+
+# Multi-error-mode helpers
+function _needs_timeseries(error_estimates)
+    any(e -> e ∈ TIMESERIES_ERRORS, error_estimates)
+end
+
+function _needs_dense(error_estimates)
+    any(e -> e ∈ DENSE_ERRORS, error_estimates)
 end
 
+"""
+    available_errors(wp_set::WorkPrecisionSet)
+
+Return the list of error estimates that were computed for this WorkPrecisionSet.
+"""
+available_errors(wp_set::WorkPrecisionSet) = wp_set.active_error_estimates
+
 function WorkPrecision(prob, alg, abstols, reltols, dts = nothing;
         name = nothing, appxsol = nothing, error_estimate = :final,
-        numruns = 20, seconds = 2, kwargs...)
+        numruns = 20, seconds = 2, timeout = nothing,
+        timeseries_errors::Union{Bool, Nothing} = nothing,
+        dense_errors::Union{Bool, Nothing} = nothing,
+        tags::Vector{Symbol} = Symbol[], kwargs...)
     N = length(abstols)
     errors = Vector{Dict{Symbol, Float64}}(undef, N)
     times = Vector{Float64}(undef, N)
@@ -201,18 +240,29 @@ function WorkPrecision(prob, alg, abstols, reltols, dts = nothing;
     end
 
     let _prob = _prob
-        timeseries_errors = error_estimate ∈ TIMESERIES_ERRORS
-        dense_errors = error_estimate ∈ DENSE_ERRORS
+        _timeseries_errors = timeseries_errors !== nothing ? timeseries_errors : (error_estimate ∈ TIMESERIES_ERRORS)
+        _dense_errors = dense_errors !== nothing ? dense_errors : (error_estimate ∈ DENSE_ERRORS)
         for i in 1:N
+            t_start = time()
             if dts === nothing
                 sol = solve(_prob, alg; kwargs..., abstol = abstols[i],
-                    reltol = reltols[i], timeseries_errors = timeseries_errors,
-                    dense_errors = dense_errors)
+                    reltol = reltols[i], timeseries_errors = _timeseries_errors,
+                    dense_errors = _dense_errors)
             else
                 sol = solve(_prob, alg; kwargs..., abstol = abstols[i],
                     reltol = reltols[i], dt = dts[i],
-                    timeseries_errors = timeseries_errors,
-                    dense_errors = dense_errors)
+                    timeseries_errors = _timeseries_errors,
+                    dense_errors = _dense_errors)
+            end
+            solve_time = time() - t_start
+
+            if timeout !== nothing && solve_time > timeout
+                @warn "Timeout exceeded for tolerance index $i ($(round(solve_time, digits=1))s > $(timeout)s)"
+                errors[i] = Dict(
+                    :l∞ => NaN, :L2 => NaN, :final => NaN, :l2 => NaN, :L∞ => NaN)
+                times[i] = NaN
+                stats[i] = nothing
+                continue
             end
 
             stats[i] = sol.stats
@@ -290,14 +340,17 @@ function WorkPrecision(prob, alg, abstols, reltols, dts = nothing;
             end
         end
     end
-    return WorkPrecision(prob, abstols, reltols, StructArray(NamedTuple.(errors)),
-        times, dts, stats, name, error_estimate, N)
+    return WorkPrecision(prob, abstols, reltols, _dicts_to_structarray(errors),
+        times, dts, stats, name, error_estimate, N, tags)
 end
 
 # Work precision information for a BVP
 function WorkPrecision(prob::AbstractBVProblem, alg, abstols, reltols, dts = nothing;
         name = nothing, appxsol = nothing, error_estimate = :final,
-        numruns = 20, seconds = 2, kwargs...)
+        numruns = 20, seconds = 2, timeout = nothing,
+        timeseries_errors::Union{Bool, Nothing} = nothing,
+        dense_errors::Union{Bool, Nothing} = nothing,
+        tags::Vector{Symbol} = Symbol[], kwargs...)
     N = length(abstols)
     errors = Vector{Dict{Symbol, Float64}}(undef, N)
     times = Vector{Float64}(undef, N)
@@ -315,18 +368,29 @@ function WorkPrecision(prob::AbstractBVProblem, alg, abstols, reltols, dts = not
     end
 
     let _prob = _prob
-        timeseries_errors = error_estimate ∈ TIMESERIES_ERRORS
-        dense_errors = error_estimate ∈ DENSE_ERRORS
+        _timeseries_errors = timeseries_errors !== nothing ? timeseries_errors : (error_estimate ∈ TIMESERIES_ERRORS)
+        _dense_errors = dense_errors !== nothing ? dense_errors : (error_estimate ∈ DENSE_ERRORS)
         for i in 1:N
+            t_start = time()
             if dts === nothing
                 sol = solve(_prob, alg; kwargs..., abstol = abstols[i],
-                    reltol = reltols[i], timeseries_errors = timeseries_errors,
-                    dense_errors = dense_errors)
+                    reltol = reltols[i], timeseries_errors = _timeseries_errors,
+                    dense_errors = _dense_errors)
             else
                 sol = solve(_prob, alg; kwargs..., abstol = abstols[i],
                     reltol = reltols[i], dt = dts[i],
-                    timeseries_errors = timeseries_errors,
-                    dense_errors = dense_errors)
+                    timeseries_errors = _timeseries_errors,
+                    dense_errors = _dense_errors)
+            end
+            solve_time = time() - t_start
+
+            if timeout !== nothing && solve_time > timeout
+                @warn "Timeout exceeded for tolerance index $i ($(round(solve_time, digits=1))s > $(timeout)s)"
+                errors[i] = Dict(
+                    :l∞ => NaN, :L2 => NaN, :final => NaN, :l2 => NaN, :L∞ => NaN)
+                times[i] = NaN
+                stats[i] = nothing
+                continue
             end
 
             stats[i] = sol.stats
@@ -347,8 +411,8 @@ function WorkPrecision(prob::AbstractBVProblem, alg, abstols, reltols, dts = not
                     end
                 else
                     errors[i] = Dict{Symbol, Float64}()
-                    for err in keys(errsol.errors)
-                        errors[i][err] = mean(errsol.errors[err])
+                    for err in keys(sol.errors)
+                        errors[i][err] = mean(sol.errors[err])
                     end
                 end
 
@@ -403,14 +467,15 @@ function WorkPrecision(prob::AbstractBVProblem, alg, abstols, reltols, dts = not
             end
         end
     end
-    return WorkPrecision(prob, abstols, reltols, StructArray(NamedTuple.(errors)),
-        times, dts, stats, name, error_estimate, N)
+    return WorkPrecision(prob, abstols, reltols, _dicts_to_structarray(errors),
+        times, dts, stats, name, error_estimate, N, tags)
 end
 
 # Work precision information for a nonlinear problem.
 function WorkPrecision(
         prob::NonlinearProblem, alg, abstols, reltols, dts = nothing; name = nothing,
-        appxsol = nothing, error_estimate = :l2, numruns = 20, seconds = 2, kwargs...)
+        appxsol = nothing, error_estimate = :l2, numruns = 20, seconds = 2,
+        timeout = nothing, tags::Vector{Symbol} = Symbol[], kwargs...)
     N = length(abstols)
     errors = Vector{Dict{Symbol, Float64}}(undef, N)
     times = Vector{Float64}(undef, N)
@@ -429,7 +494,17 @@ function WorkPrecision(
 
     let _prob = _prob
         for i in 1:N
+            t_start = time()
             sol = solve(_prob, alg; kwargs..., abstol = abstols[i], reltol = reltols[i])
+            solve_time = time() - t_start
+
+            if timeout !== nothing && solve_time > timeout
+                @warn "Timeout exceeded for tolerance index $i ($(round(solve_time, digits=1))s > $(timeout)s)"
+                errors[i] = Dict(error_estimate => NaN)
+                times[i] = NaN
+                stats[i] = nothing
+                continue
+            end
 
             stats[i] = sol.stats
 
@@ -461,35 +536,44 @@ function WorkPrecision(
         end
     end
 
-    return WorkPrecision(prob, abstols, reltols, StructArray(NamedTuple.(errors)),
-        times, dts, stats, name, error_estimate, N)
+    return WorkPrecision(prob, abstols, reltols, _dicts_to_structarray(errors),
+        times, dts, stats, name, error_estimate, N, tags)
 end
 
 function WorkPrecisionSet(prob,
         abstols, reltols, setups;
         print_names = false, names = nothing, appxsol = nothing,
-        error_estimate = :final,
-        test_dt = nothing, kwargs...)
+        error_estimate = :final, error_estimates = nothing,
+        test_dt = nothing, timeout = nothing, kwargs...)
     N = length(setups)
     @assert names === nothing || length(setups) == length(names)
     wps = Vector{WorkPrecision}(undef, N)
     if names === nothing
         names = [_default_name(setup[:alg]) for setup in setups]
     end
+
+    _active = error_estimates !== nothing ? collect(error_estimates) : [error_estimate]
+    _ts_errors = error_estimates !== nothing ? _needs_timeseries(_active) : nothing
+    _dense_errors = error_estimates !== nothing ? _needs_dense(_active) : nothing
+
     for i in 1:N
         print_names && println(names[i])
         _abstols = get(setups[i], :abstols, abstols)
         _reltols = get(setups[i], :reltols, reltols)
         _dts = get(setups[i], :dts, nothing)
+        _tags = get(setups[i], :tags, Symbol[])
         filtered_setup = filter(p -> p.first in DiffEqBase.allowedkeywords, setups[i])
 
         wps[i] = WorkPrecision(prob, setups[i][:alg], _abstols, _reltols, _dts;
             appxsol = appxsol,
             error_estimate = error_estimate,
-            name = names[i], kwargs..., filtered_setup...)
+            name = names[i], timeout = timeout,
+            timeseries_errors = _ts_errors,
+            dense_errors = _dense_errors,
+            tags = _tags, kwargs..., filtered_setup...)
     end
     return WorkPrecisionSet(wps, N, abstols, reltols, prob, setups, names, error_estimate,
-        nothing)
+        nothing, _active)
 end
 
 @def error_calculation begin
@@ -540,7 +624,8 @@ function WorkPrecisionSet(prob::AbstractRODEProblem, abstols, reltols, setups,
         test_dt = nothing;
         numruns = 20, numruns_error = 20,
         print_names = false, names = nothing, appxsol_setup = nothing,
-        error_estimate = :final, parallel_type = :none,
+        error_estimate = :final, error_estimates = nothing,
+        parallel_type = :none, timeout = nothing,
         kwargs...)
     @assert names === nothing || length(setups) == length(names)
     timeseries_errors = DiffEqBase.has_analytic(prob.f) &&
@@ -557,6 +642,8 @@ function WorkPrecisionSet(prob::AbstractRODEProblem, abstols, reltols, setups,
     end
     time_tmp = Vector{Float64}(undef, numruns)
 
+    _active = error_estimates !== nothing ? collect(error_estimates) : [error_estimate]
+
     # First calculate all of the errors
     if parallel_type == :threads
         Threads.@threads for i in 1:numruns_error
@@ -600,26 +687,34 @@ function WorkPrecisionSet(prob::AbstractRODEProblem, abstols, reltols, setups,
         x = isempty(_sol.t) ? 0 : round(Int, mean(_sol.t) - sum(_sol.t) / length(_sol.t))
         GC.gc()
         for j in 1:M
+            timed_out = false
             for i in 1:numruns
+                t_start = time()
                 time_tmp[i] = @elapsed sol = solve(prob, setups[k][:alg];
                     kwargs..., filtered_setup...,
                     abstol = _abstols[k][j],
                     reltol = _reltols[k][j], dt = _dts[k][j],
                     timeseries_errors = false,
                     dense_errors = false)
+                if timeout !== nothing && (time() - t_start) > timeout
+                    @warn "Timeout exceeded for method $k, tolerance $j"
+                    timed_out = true
+                    break
+                end
             end
-            times[j, k] = mean(time_tmp) + x
+            times[j, k] = timed_out ? NaN : mean(time_tmp) + x
             GC.gc()
         end
     end
 
     stats = nothing
     wps = [WorkPrecision(prob, _abstols[i], _reltols[i],
-               StructArray(NamedTuple.(errors[i])),
-               times[:, i], _dts[i], stats, names[i], error_estimate, N)
+               _dicts_to_structarray(errors[i]),
+               times[:, i], _dts[i], stats, names[i], error_estimate, N,
+               get(setups[i], :tags, Symbol[]))
            for i in 1:N]
     WorkPrecisionSet(wps, N, abstols, reltols, prob, setups, names, error_estimate,
-        numruns_error)
+        numruns_error, _active)
 end
 
 function WorkPrecisionSet(prob::AbstractEnsembleProblem, abstols, reltols, setups,
@@ -627,8 +722,9 @@ function WorkPrecisionSet(prob::AbstractEnsembleProblem, abstols, reltols, setup
         numruns = 5, trajectories = 1000,
         print_names = false, names = nothing, appxsol_setup = nothing,
         expected_value = nothing,
-        error_estimate = :weak_final, ensemblealg = EnsembleThreads(),
-        kwargs...)
+        error_estimate = :weak_final, error_estimates = nothing,
+        ensemblealg = EnsembleThreads(),
+        timeout = nothing, kwargs...)
     @assert names === nothing || length(setups) == length(names)
 
     weak_timeseries_errors = error_estimate ∈ WEAK_TIMESERIES_ERRORS
@@ -643,6 +739,8 @@ function WorkPrecisionSet(prob::AbstractEnsembleProblem, abstols, reltols, setup
     end
     time_tmp = Vector{Float64}(undef, numruns)
 
+    _active = error_estimates !== nothing ? collect(error_estimates) : [error_estimate]
+
     # First calculate all of the errors
     _abstols = [get(setups[k], :abstols, abstols) for k in 1:N]
     _reltols = [get(setups[k], :reltols, reltols) for k in 1:N]
@@ -708,7 +806,9 @@ function WorkPrecisionSet(prob::AbstractEnsembleProblem, abstols, reltols, setup
         #x = isempty(_sol.t) ? 0 : round(Int,mean(_sol.t) - sum(_sol.t)/length(_sol.t))
         GC.gc()
         for j in 1:M
+            timed_out = false
             for i in 1:numruns
+                t_start = time()
                 time_tmp[i] = @elapsed sol = solve(prob, setups[k][:alg], ensemblealg;
                     filtered_setup...,
                     abstol = _abstols[k][j],
@@ -718,44 +818,59 @@ function WorkPrecisionSet(prob::AbstractEnsembleProblem, abstols, reltols, setup
                     dense_errors = false,
                     trajectories = Int(trajectories),
                     kwargs...)
+                if timeout !== nothing && (time() - t_start) > timeout
+                    @warn "Timeout exceeded for method $k, tolerance $j"
+                    timed_out = true
+                    break
+                end
             end
-            times[j, k] = mean(time_tmp) #+ x
+            times[j, k] = timed_out ? NaN : mean(time_tmp) #+ x
             GC.gc()
         end
     end
     stats = nothing
     wps = [WorkPrecision(prob, _abstols[i], _reltols[i], errors[i], times[:, i],
-               _dts[i], stats, names[i], error_estimate, N)
+               _dts[i], stats, names[i], error_estimate, N,
+               get(setups[i], :tags, Symbol[]))
            for i in 1:N]
     WorkPrecisionSet(wps, N, abstols, reltols, prob, setups, names, error_estimate,
-        Int(trajectories))
+        Int(trajectories), _active)
 end
 
 function WorkPrecisionSet(prob::AbstractBVProblem,
         abstols, reltols, setups;
         print_names = false, names = nothing, appxsol = nothing,
-        error_estimate = :final,
-        test_dt = nothing, kwargs...)
+        error_estimate = :final, error_estimates = nothing,
+        test_dt = nothing, timeout = nothing, kwargs...)
     N = length(setups)
     @assert names === nothing || length(setups) == length(names)
     wps = Vector{WorkPrecision}(undef, N)
     if names === nothing
         names = [_default_name(setup[:alg]) for setup in setups]
     end
+
+    _active = error_estimates !== nothing ? collect(error_estimates) : [error_estimate]
+    _ts_errors = error_estimates !== nothing ? _needs_timeseries(_active) : nothing
+    _dense_errors = error_estimates !== nothing ? _needs_dense(_active) : nothing
+
     for i in 1:N
         print_names && println(names[i])
         _abstols = get(setups[i], :abstols, abstols)
         _reltols = get(setups[i], :reltols, reltols)
         _dts = get(setups[i], :dts, nothing)
+        _tags = get(setups[i], :tags, Symbol[])
         filtered_setup = filter(p -> p.first in DiffEqBase.allowedkeywords, setups[i])
 
         wps[i] = WorkPrecision(prob, setups[i][:alg], _abstols, _reltols, _dts;
             appxsol = appxsol,
             error_estimate = error_estimate,
-            name = names[i], kwargs..., filtered_setup...)
+            name = names[i], timeout = timeout,
+            timeseries_errors = _ts_errors,
+            dense_errors = _dense_errors,
+            tags = _tags, kwargs..., filtered_setup...)
     end
     return WorkPrecisionSet(wps, N, abstols, reltols, prob, setups, names, error_estimate,
-        nothing)
+        nothing, _active)
 end
 
 function get_sample_errors(prob::AbstractRODEProblem, setup, test_dt = nothing;
@@ -823,6 +938,185 @@ function get_sample_errors(prob::AbstractRODEProblem, setup, test_dt = nothing;
     end
 end
 
+## Tagging and filtering helpers
+
+"""
+    filter_by_tags(wp_set::WorkPrecisionSet, tags::Symbol...) -> WorkPrecisionSet
+
+Return a new `WorkPrecisionSet` containing only entries whose tags include
+ALL of the specified tags (AND logic).
+"""
+function filter_by_tags(wp_set::WorkPrecisionSet, tags::Symbol...)
+    isempty(tags) && return wp_set
+    indices = findall(wp -> all(t -> t in wp.tags, tags), wp_set.wps)
+    _subset_wps(wp_set, indices)
+end
+
+"""
+    exclude_by_tags(wp_set::WorkPrecisionSet, tags::Symbol...) -> WorkPrecisionSet
+
+Return a new `WorkPrecisionSet` excluding entries that have ANY of the specified tags.
+"""
+function exclude_by_tags(wp_set::WorkPrecisionSet, tags::Symbol...)
+    isempty(tags) && return wp_set
+    indices = findall(wp -> !any(t -> t in wp.tags, tags), wp_set.wps)
+    _subset_wps(wp_set, indices)
+end
+
+"""
+    get_tags(wp_set::WorkPrecisionSet) -> Vector{Vector{Symbol}}
+
+Return the tags for each entry in the `WorkPrecisionSet`.
+"""
+get_tags(wp_set::WorkPrecisionSet) = [wp.tags for wp in wp_set.wps]
+
+"""
+    unique_tags(wp_set::WorkPrecisionSet) -> Vector{Symbol}
+
+Return all unique tags present across entries in the `WorkPrecisionSet`.
+"""
+function unique_tags(wp_set::WorkPrecisionSet)
+    alltags = Symbol[]
+    for wp in wp_set.wps
+        append!(alltags, wp.tags)
+    end
+    return unique!(sort!(alltags))
+end
+
+"""
+    merge_wp_sets(sets::WorkPrecisionSet...) -> WorkPrecisionSet
+
+Merge multiple `WorkPrecisionSet`s into a single set. All entries are combined.
+The metadata (abstols, reltols, prob, error_estimate) is taken from the first set.
+"""
+function merge_wp_sets(sets::WorkPrecisionSet...)
+    isempty(sets) && throw(ArgumentError("At least one WorkPrecisionSet is required"))
+    wps = vcat([s.wps for s in sets]...)
+    all_setups = vcat([s.setups for s in sets]...)
+    all_names = vcat([s.names for s in sets]...)
+    N = length(wps)
+    first_set = first(sets)
+    return WorkPrecisionSet(
+        wps, N, first_set.abstols, first_set.reltols, first_set.prob,
+        all_setups, all_names, first_set.error_estimate, first_set.numruns)
+end
+
+function _subset_wps(wp_set::WorkPrecisionSet, indices::Vector{Int})
+    isempty(indices) &&
+        @warn "No entries match the specified tags. Returning empty WorkPrecisionSet."
+    wps = wp_set.wps[indices]
+    setups = wp_set.setups isa AbstractVector ? wp_set.setups[indices] : wp_set.setups
+    names = wp_set.names isa AbstractVector ? wp_set.names[indices] : wp_set.names
+    return WorkPrecisionSet(
+        wps, length(wps), wp_set.abstols, wp_set.reltols, wp_set.prob,
+        setups, names, wp_set.error_estimate, wp_set.numruns)
+end
+
+## Best-of-family helpers
+
+"""
+    wp_area(wp::WorkPrecision)
+
+Compute the area under the log-log work-precision curve using trapezoidal integration.
+Lower area = better (less time for given error). Returns `Inf` if fewer than 2 valid points.
+"""
+function wp_area(wp::WorkPrecision)
+    errs = getproperty(wp.errors, wp.error_estimate)
+    times = wp.times
+    valid = [(e, t) for (e, t) in zip(errs, times) if !isnan(e) && !isnan(t) && e > 0 && t > 0]
+    length(valid) < 2 && return Inf
+    log_errs = [log10(v[1]) for v in valid]
+    log_times = [log10(v[2]) for v in valid]
+    perm = sortperm(log_errs)
+    log_errs = log_errs[perm]
+    log_times = log_times[perm]
+    area = 0.0
+    for i in 2:length(log_errs)
+        area += 0.5 * (log_times[i] + log_times[i - 1]) * (log_errs[i] - log_errs[i - 1])
+    end
+    return area
+end
+
+"""
+    best_by_tag(wp_set, tag; n=1, metric=:area)
+
+Return the top `n` methods matching `tag`, ranked by work-precision performance.
+"""
+function best_by_tag(wp_set::WorkPrecisionSet, tag::Symbol; n::Int = 1, metric::Symbol = :area)
+    filtered = filter_by_tags(wp_set, tag)
+    length(filtered) == 0 && return filtered
+    if metric == :area
+        areas = [wp_area(wp) for wp in filtered.wps]
+        perm = sortperm(areas)
+        selected = perm[1:min(n, length(perm))]
+        return _subset_wps(filtered, selected)
+    else
+        throw(ArgumentError("Unknown metric: $metric. Supported: :area"))
+    end
+end
+
+"""
+    best_of_families(wp_set, family_tags; n=1, metric=:area)
+
+Select the best `n` methods from each family tag and combine into a single WorkPrecisionSet.
+"""
+function best_of_families(wp_set::WorkPrecisionSet, family_tags; n::Int = 1, metric::Symbol = :area)
+    results = WorkPrecisionSet[]
+    for tag in family_tags
+        best = best_by_tag(wp_set, tag; n = n, metric = metric)
+        length(best) > 0 && push!(results, best)
+    end
+    isempty(results) && throw(ArgumentError("No methods found for any of the specified family tags"))
+    return merge_wp_sets(results...)
+end
+
+## AutoDiff comparison helpers
+
+"""
+    ad_backend_name(backend)
+
+Get a short string name from an AD backend object. Strips "Auto" prefix.
+"""
+function ad_backend_name(backend)
+    name = string(nameof(typeof(backend)))
+    startswith(name, "Auto") ? name[5:end] : name
+end
+
+"""
+    with_autodiff_variants(setups; ad_backends, tag_prefix=:autodiff)
+
+Create AD variants of setup dicts. For each setup and each AD backend, creates a
+new setup with the `:autodiff` key set and tags augmented with an AD-specific tag.
+Original setups get tagged with `Symbol(tag_prefix, "_default")`.
+"""
+function with_autodiff_variants(setups; ad_backends, tag_prefix::Symbol = :autodiff)
+    result = Vector{Dict{Symbol, Any}}()
+
+    for setup in setups
+        # Tag original setup (copy to avoid mutation)
+        original = copy(setup)
+        original_tags = copy(get(original, :tags, Symbol[]))
+        push!(original_tags, Symbol(tag_prefix, :_default))
+        original[:tags] = original_tags
+        push!(result, original)
+
+        # Create variants for each AD backend
+        for backend in ad_backends
+            variant = copy(setup)
+            variant[:autodiff] = backend
+            variant_tags = copy(get(variant, :tags, Symbol[]))
+            backend_name = lowercase(ad_backend_name(backend))
+            push!(variant_tags, Symbol(tag_prefix, :_, Symbol(backend_name)))
+            variant[:tags] = variant_tags
+            push!(result, variant)
+        end
+    end
+
+    return result
+end
+
+## Base overloads
+
 Base.length(wp::WorkPrecision) = wp.N
 Base.size(wp::WorkPrecision) = length(wp)
 Base.getindex(wp::WorkPrecision, i::Int) = wp.times[i]
diff --git a/lib/DiffEqDevTools/src/plotrecipes.jl b/lib/DiffEqDevTools/src/plotrecipes.jl
index e0018800a12..54aa66f17c3 100644
--- a/lib/DiffEqDevTools/src/plotrecipes.jl
+++ b/lib/DiffEqDevTools/src/plotrecipes.jl
@@ -91,7 +91,7 @@ function key_to_label(key::Symbol)
     elseif key == :maxeig
         return "Maximum eigenvalue recorded"
     else
-        return key
+        return String(key)
     end
 end
 
@@ -99,20 +99,90 @@ end
         x = wp_set.error_estimate,
         y = :times,
         view = :benchmark,
-        color = nothing)
+        color = nothing,
+        tags = nothing,
+        include_tags = nothing,
+        exclude_tags = nothing,
+        reference_tags = nothing,
+        reference_style = (linestyle = :dash, linewidth = 1, alpha = 0.5))
+    # Apply tag-based filtering if specified
+    if tags !== nothing || include_tags !== nothing || exclude_tags !== nothing
+        filtered = wp_set
+        if tags !== nothing
+            filtered = filter_by_tags(filtered, tags...)
+        end
+        if include_tags !== nothing
+            extra = filter_by_tags(wp_set, include_tags...)
+            seen = Set(wp.name for wp in filtered.wps)
+            extra_wps = [wp for wp in extra.wps if wp.name ∉ seen]
+            if !isempty(extra_wps)
+                all_wps = vcat(filtered.wps, extra_wps)
+                all_names = [wp.name for wp in all_wps]
+                filtered = WorkPrecisionSet(
+                    all_wps, length(all_wps), filtered.abstols, filtered.reltols,
+                    filtered.prob, filtered.setups, all_names, filtered.error_estimate,
+                    filtered.numruns)
+            end
+        end
+        if exclude_tags !== nothing
+            filtered = exclude_by_tags(filtered, exclude_tags...)
+        end
+        wp_set = filtered
+    end
+
     if view == :benchmark
-        seriestype --> :path
-        linewidth --> 3
-        xscale --> :log10
-        yscale --> :log10
-        markershape --> :auto
-        xs = [get_val_from_wp(wp, x) for wp in wp_set.wps]
-        ys = [get_val_from_wp(wp, y) for wp in wp_set.wps]
-        xguide --> key_to_label(x)
-        yguide --> key_to_label(y)
-        legend --> :outerright
-        label --> reshape(wp_set.names, 1, length(wp_set))
-        return xs, ys
+        if reference_tags !== nothing
+            ref_tags = reference_tags isa Symbol ? (reference_tags,) : reference_tags
+            ref_indices = findall(wp -> any(t -> t in wp.tags, ref_tags), wp_set.wps)
+            main_indices = setdiff(1:length(wp_set.wps), ref_indices)
+
+            if !isempty(ref_indices)
+                @series begin
+                    seriestype --> :path
+                    linewidth --> get(reference_style, :linewidth, 1)
+                    linestyle --> get(reference_style, :linestyle, :dash)
+                    seriesalpha --> get(reference_style, :alpha, 0.5)
+                    xscale --> :log10
+                    yscale --> :log10
+                    markershape --> :auto
+                    xguide --> key_to_label(x)
+                    yguide --> key_to_label(y)
+                    legend --> :outerright
+                    ref_xs = [get_val_from_wp(wp_set.wps[i], x) for i in ref_indices]
+                    ref_ys = [get_val_from_wp(wp_set.wps[i], y) for i in ref_indices]
+                    label --> reshape([wp_set.wps[i].name for i in ref_indices], 1, length(ref_indices))
+                    ref_xs, ref_ys
+                end
+            end
+
+            if !isempty(main_indices)
+                seriestype --> :path
+                linewidth --> 3
+                xscale --> :log10
+                yscale --> :log10
+                markershape --> :auto
+                xguide --> key_to_label(x)
+                yguide --> key_to_label(y)
+                legend --> :outerright
+                main_xs = [get_val_from_wp(wp_set.wps[i], x) for i in main_indices]
+                main_ys = [get_val_from_wp(wp_set.wps[i], y) for i in main_indices]
+                label --> reshape([wp_set.wps[i].name for i in main_indices], 1, length(main_indices))
+                return main_xs, main_ys
+            end
+        else
+            seriestype --> :path
+            linewidth --> 3
+            xscale --> :log10
+            yscale --> :log10
+            markershape --> :auto
+            xs = [get_val_from_wp(wp, x) for wp in wp_set.wps]
+            ys = [get_val_from_wp(wp, y) for wp in wp_set.wps]
+            xguide --> key_to_label(x)
+            yguide --> key_to_label(y)
+            legend --> :outerright
+            label --> reshape(wp_set.names, 1, length(wp_set))
+            return xs, ys
+        end
     elseif view == :dt_convergence
         idts = filter(i -> haskey(wp_set.setups[i], :dts), 1:length(wp_set))
         length(idts) > 0 ||
diff --git a/lib/DiffEqDevTools/test/autodiff_tests.jl b/lib/DiffEqDevTools/test/autodiff_tests.jl
new file mode 100644
index 00000000000..ab98b1dbb00
--- /dev/null
+++ b/lib/DiffEqDevTools/test/autodiff_tests.jl
@@ -0,0 +1,54 @@
+using DiffEqDevTools, Test
+
+# Mock AD backend types for testing
+struct AutoForwardDiff end
+struct AutoFiniteDiff end
+struct ManualJac end
+
+@testset "AutoDiff Helpers" begin
+    @testset "ad_backend_name" begin
+        @test DiffEqDevTools.ad_backend_name(AutoForwardDiff()) == "ForwardDiff"
+        @test DiffEqDevTools.ad_backend_name(AutoFiniteDiff()) == "FiniteDiff"
+        @test DiffEqDevTools.ad_backend_name(ManualJac()) == "ManualJac"
+    end
+
+    @testset "with_autodiff_variants count" begin
+        setups = [
+            Dict{Symbol, Any}(:alg => nothing, :tags => [:rk]),
+            Dict{Symbol, Any}(:alg => nothing),
+        ]
+        expanded = with_autodiff_variants(setups;
+            ad_backends = [AutoForwardDiff(), AutoFiniteDiff()])
+        # 2 originals + 2x2 variants = 6
+        @test length(expanded) == 6
+    end
+
+    @testset "original gets default tag" begin
+        setups = [Dict{Symbol, Any}(:alg => nothing, :tags => [:rk])]
+        expanded = with_autodiff_variants(setups; ad_backends = [AutoForwardDiff()])
+        @test :autodiff_default in expanded[1][:tags]
+        @test :rk in expanded[1][:tags]
+    end
+
+    @testset "variant gets backend tag and autodiff key" begin
+        setups = [Dict{Symbol, Any}(:alg => nothing, :tags => [:rk])]
+        expanded = with_autodiff_variants(setups; ad_backends = [AutoForwardDiff()])
+        @test expanded[2][:autodiff] isa AutoForwardDiff
+        @test :rk in expanded[2][:tags]
+    end
+
+    @testset "custom tag_prefix" begin
+        setups = [Dict{Symbol, Any}(:alg => nothing)]
+        expanded = with_autodiff_variants(setups;
+            ad_backends = [AutoForwardDiff()], tag_prefix = :ad)
+        @test :ad_default in expanded[1][:tags]
+    end
+
+    @testset "no mutation of original" begin
+        original_setup = Dict{Symbol, Any}(:alg => nothing, :tags => [:rk])
+        setups = [original_setup]
+        expanded = with_autodiff_variants(setups; ad_backends = [AutoForwardDiff()])
+        @test :autodiff_default ∉ original_setup[:tags]
+        @test length(original_setup[:tags]) == 1
+    end
+end
diff --git a/lib/DiffEqDevTools/test/autoplot_tests.jl b/lib/DiffEqDevTools/test/autoplot_tests.jl
new file mode 100644
index 00000000000..f9fdf677e20
--- /dev/null
+++ b/lib/DiffEqDevTools/test/autoplot_tests.jl
@@ -0,0 +1,49 @@
+using OrdinaryDiffEq, DiffEqDevTools, Test
+using ODEProblemLibrary: prob_ode_linear
+
+prob = prob_ode_linear
+
+@testset "Autoplot" begin
+    abstols = 1 ./ 10 .^ (3:6)
+    reltols = 1 ./ 10 .^ (3:6)
+
+    setups = [
+        Dict{Symbol, Any}(:alg => RK4(), :tags => [:rk, :fourth_order]),
+        Dict{Symbol, Any}(:alg => DP5(), :tags => [:rk, :fifth_order]),
+        Dict{Symbol, Any}(:alg => Tsit5(), :tags => [:rk, :fifth_order]),
+        Dict{Symbol, Any}(:alg => BS3(), :tags => [:rk, :third_order]),
+        Dict{Symbol, Any}(:alg => Euler(), :tags => [:euler, :first_order, :reference]),
+    ]
+
+    wp_set = WorkPrecisionSet(prob, abstols, reltols, setups; dt = 1 / 2^4, numruns = 2)
+
+    @testset "explicit families" begin
+        plots = autoplot(wp_set; families = [:rk, :euler])
+        @test haskey(plots, "family_rk")
+        @test haskey(plots, "family_euler")
+        @test haskey(plots, "best_of_families")
+        @test haskey(plots, "all")
+        @test length(plots["family_rk"]) == 4  # RK4, DP5, Tsit5, BS3
+        @test length(plots["family_euler"]) == 1
+    end
+
+    @testset "auto-detect families" begin
+        plots = autoplot(wp_set)
+        @test haskey(plots, "all")
+        # Should detect families that aren't on >80% of methods
+    end
+
+    @testset "with reference_tags" begin
+        plots = autoplot(wp_set; families = [:rk], reference_tags = :reference)
+        # rk family should include the reference method too
+        rk_names = plots["family_rk"].names
+        @test "Euler" in rk_names || length(plots["family_rk"]) >= 4
+    end
+
+    @testset "best_n" begin
+        plots = autoplot(wp_set; families = [:rk, :euler], best_n = 1)
+        bof = plots["best_of_families"]
+        # Should have at most 1 from each family
+        @test length(bof) <= 2
+    end
+end
diff --git a/lib/DiffEqDevTools/test/best_of_family_tests.jl b/lib/DiffEqDevTools/test/best_of_family_tests.jl
new file mode 100644
index 00000000000..a8c405f438e
--- /dev/null
+++ b/lib/DiffEqDevTools/test/best_of_family_tests.jl
@@ -0,0 +1,50 @@
+using OrdinaryDiffEq, DiffEqDevTools, Test
+using ODEProblemLibrary: prob_ode_linear
+
+prob = prob_ode_linear
+
+@testset "Best-of-Family" begin
+    abstols = 1 ./ 10 .^ (3:6)
+    reltols = 1 ./ 10 .^ (3:6)
+
+    setups = [
+        Dict{Symbol, Any}(:alg => RK4(), :tags => [:rk, :fourth_order]),
+        Dict{Symbol, Any}(:alg => DP5(), :tags => [:rk, :fifth_order]),
+        Dict{Symbol, Any}(:alg => Tsit5(), :tags => [:rk, :fifth_order]),
+        Dict{Symbol, Any}(:alg => BS3(), :tags => [:rk, :third_order]),
+        Dict{Symbol, Any}(:alg => Euler(), :tags => [:euler, :first_order]),
+    ]
+
+    wp_set = WorkPrecisionSet(prob, abstols, reltols, setups; dt = 1 / 2^4, numruns = 2)
+
+    @testset "wp_area" begin
+        area = wp_area(wp_set[1])
+        @test isfinite(area)
+    end
+
+    @testset "best_by_tag" begin
+        best_rk = best_by_tag(wp_set, :rk; n = 2)
+        @test length(best_rk) == 2
+
+        best_euler = best_by_tag(wp_set, :euler; n = 1)
+        @test length(best_euler) == 1
+        @test best_euler.names == ["Euler"]
+
+        # No match
+        empty = best_by_tag(wp_set, :nonexistent)
+        @test length(empty) == 0
+    end
+
+    @testset "best_of_families" begin
+        best = best_of_families(wp_set, [:rk, :euler]; n = 1)
+        @test length(best) == 2
+        @test "Euler" in best.names
+
+        # Error on all-empty
+        @test_throws ArgumentError best_of_families(wp_set, [:nonexistent])
+    end
+
+    @testset "unknown metric" begin
+        @test_throws ArgumentError best_by_tag(wp_set, :rk; metric = :unknown)
+    end
+end
diff --git a/lib/DiffEqDevTools/test/multi_error_tests.jl b/lib/DiffEqDevTools/test/multi_error_tests.jl
new file mode 100644
index 00000000000..ca9aea07b28
--- /dev/null
+++ b/lib/DiffEqDevTools/test/multi_error_tests.jl
@@ -0,0 +1,31 @@
+using OrdinaryDiffEq, DiffEqDevTools, Test
+using ODEProblemLibrary: prob_ode_linear
+
+prob = prob_ode_linear
+
+@testset "Multi-error-mode" begin
+    abstols = 1 ./ 10 .^ (3:6)
+    reltols = 1 ./ 10 .^ (3:6)
+
+    @testset "Backward compatible 9-arg constructor" begin
+        wp = WorkPrecision(prob, DP5(), abstols, reltols; name = "DP5", numruns = 2)
+        wps = WorkPrecisionSet([wp], 1, abstols, reltols, prob, [], ["DP5"], :final, nothing)
+        @test wps.active_error_estimates == [:final]
+        @test available_errors(wps) == [:final]
+    end
+
+    @testset "Single error_estimate (default)" begin
+        setups = [Dict{Symbol, Any}(:alg => DP5())]
+        wp_set = WorkPrecisionSet(prob, abstols, reltols, setups; numruns = 2)
+        @test available_errors(wp_set) == [:final]
+    end
+
+    @testset "Multi error_estimates" begin
+        setups = [Dict{Symbol, Any}(:alg => DP5())]
+        wp_set = WorkPrecisionSet(prob, abstols, reltols, setups;
+            numruns = 2, error_estimates = [:final, :l2])
+        @test Set(available_errors(wp_set)) == Set([:final, :l2])
+        # l2 errors should be computed (timeseries_errors was enabled)
+        @test :l2 in propertynames(wp_set[1].errors)
+    end
+end
diff --git a/lib/DiffEqDevTools/test/runtests.jl b/lib/DiffEqDevTools/test/runtests.jl
index d6447a41069..2c9e027be5e 100644
--- a/lib/DiffEqDevTools/test/runtests.jl
+++ b/lib/DiffEqDevTools/test/runtests.jl
@@ -26,3 +26,21 @@ end
 @time @testset "Plot Recipes (Nonlinearsolve WP-diagrams)" begin
     include("nonlinearsolve_wpdiagram_tests.jl")
 end
+@time @testset "Tagging Tests" begin
+    include("tagging_tests.jl")
+end
+@time @testset "Multi-Error Tests" begin
+    include("multi_error_tests.jl")
+end
+@time @testset "Best-of-Family Tests" begin
+    include("best_of_family_tests.jl")
+end
+@time @testset "Timeout Tests" begin
+    include("timeout_tests.jl")
+end
+@time @testset "AutoDiff Helper Tests" begin
+    include("autodiff_tests.jl")
+end
+@time @testset "Autoplot Tests" begin
+    include("autoplot_tests.jl")
+end
diff --git a/lib/DiffEqDevTools/test/tagging_tests.jl b/lib/DiffEqDevTools/test/tagging_tests.jl
new file mode 100644
index 00000000000..aa57b132c0c
--- /dev/null
+++ b/lib/DiffEqDevTools/test/tagging_tests.jl
@@ -0,0 +1,150 @@
+using OrdinaryDiffEq, DiffEqDevTools, DiffEqBase, Test
+using ODEProblemLibrary: prob_ode_linear
+
+prob = prob_ode_linear
+
+@testset "Tags in WorkPrecision" begin
+    abstols = 1 ./ 10 .^ (3:6)
+    reltols = 1 ./ 10 .^ (3:6)
+
+    tgs = Vector{Symbol}([:rk, :fifth_order])
+    wp = WorkPrecision(
+        prob, DP5(), abstols, reltols;
+        name = "DP5", tags = tgs
+    )
+    @test wp.tags == [:rk, :fifth_order]
+
+    wp_notags = WorkPrecision(prob, DP5(), abstols, reltols; name = "DP5")
+    @test isempty(wp_notags.tags)
+end
+
+@testset "Tags in WorkPrecisionSet via setups" begin
+    abstols = 1 ./ 10 .^ (3:6)
+    reltols = 1 ./ 10 .^ (3:6)
+
+    setups = [
+        Dict{Symbol, Any}(:alg => RK4(), :tags => [:rk, :fourth_order]),
+        Dict{Symbol, Any}(:alg => DP5(), :tags => [:rk, :fifth_order]),
+        Dict{Symbol, Any}(:alg => BS3(), :tags => [:rk, :third_order]),
+        Dict{Symbol, Any}(:alg => Euler(), :tags => [:rk, :first_order, :reference]),
+        Dict{Symbol, Any}(:alg => Midpoint()),
+    ]
+
+    wp_set = WorkPrecisionSet(
+        prob, abstols, reltols, setups; dt = 1 / 2^4, numruns = 2
+    )
+
+    @test wp_set[1].tags == [:rk, :fourth_order]
+    @test wp_set[2].tags == [:rk, :fifth_order]
+    @test wp_set[3].tags == [:rk, :third_order]
+    @test wp_set[4].tags == [:rk, :first_order, :reference]
+    @test isempty(wp_set[5].tags)
+end
+
+@testset "filter_by_tags" begin
+    abstols = 1 ./ 10 .^ (3:6)
+    reltols = 1 ./ 10 .^ (3:6)
+
+    setups = [
+        Dict{Symbol, Any}(:alg => RK4(), :tags => [:rk, :fourth_order]),
+        Dict{Symbol, Any}(:alg => DP5(), :tags => [:rk, :fifth_order]),
+        Dict{Symbol, Any}(:alg => BS3(), :tags => [:rk, :third_order]),
+        Dict{Symbol, Any}(:alg => Euler(), :tags => [:euler, :first_order, :reference]),
+        Dict{Symbol, Any}(:alg => Midpoint()),
+    ]
+
+    wp_set = WorkPrecisionSet(
+        prob, abstols, reltols, setups; dt = 1 / 2^4, numruns = 2
+    )
+
+    # Filter by single tag
+    rk_only = filter_by_tags(wp_set, :rk)
+    @test length(rk_only) == 3
+    @test Set(rk_only.names) == Set(["RK4", "DP5", "BS3"])
+
+    # Filter by multiple tags (AND logic)
+    rk_5th = filter_by_tags(wp_set, :rk, :fifth_order)
+    @test length(rk_5th) == 1
+    @test rk_5th.names == ["DP5"]
+
+    # Filter by tag that matches nothing
+    empty_result = filter_by_tags(wp_set, :nonexistent)
+    @test length(empty_result) == 0
+
+    # Empty tags returns original
+    same = filter_by_tags(wp_set)
+    @test length(same) == length(wp_set)
+end
+
+@testset "exclude_by_tags" begin
+    abstols = 1 ./ 10 .^ (3:6)
+    reltols = 1 ./ 10 .^ (3:6)
+
+    setups = [
+        Dict{Symbol, Any}(:alg => RK4(), :tags => [:rk, :fourth_order]),
+        Dict{Symbol, Any}(:alg => DP5(), :tags => [:rk, :fifth_order]),
+        Dict{Symbol, Any}(:alg => Euler(), :tags => [:euler, :reference]),
+        Dict{Symbol, Any}(:alg => Midpoint()),
+    ]
+
+    wp_set = WorkPrecisionSet(
+        prob, abstols, reltols, setups; dt = 1 / 2^4, numruns = 2
+    )
+
+    # Exclude reference
+    no_ref = exclude_by_tags(wp_set, :reference)
+    @test length(no_ref) == 3
+    @test "Euler" ∉ no_ref.names
+
+    # Exclude multiple (OR logic)
+    no_rk_or_ref = exclude_by_tags(wp_set, :rk, :reference)
+    @test length(no_rk_or_ref) == 1
+    @test no_rk_or_ref.names == ["Midpoint"]
+
+    # Empty exclude returns original
+    same = exclude_by_tags(wp_set)
+    @test length(same) == length(wp_set)
+end
+
+@testset "unique_tags and get_tags" begin
+    abstols = 1 ./ 10 .^ (3:6)
+    reltols = 1 ./ 10 .^ (3:6)
+
+    setups = [
+        Dict{Symbol, Any}(:alg => RK4(), :tags => [:rk, :fourth_order]),
+        Dict{Symbol, Any}(:alg => DP5(), :tags => [:rk, :fifth_order]),
+        Dict{Symbol, Any}(:alg => Euler(), :tags => [:euler, :reference]),
+    ]
+
+    wp_set = WorkPrecisionSet(
+        prob, abstols, reltols, setups; dt = 1 / 2^4, numruns = 2
+    )
+
+    tags = get_tags(wp_set)
+    @test tags[1] == [:rk, :fourth_order]
+    @test tags[2] == [:rk, :fifth_order]
+    @test tags[3] == [:euler, :reference]
+
+    utags = unique_tags(wp_set)
+    @test Set(utags) == Set([:rk, :fourth_order, :fifth_order, :euler, :reference])
+end
+
+@testset "merge_wp_sets" begin
+    abstols = 1 ./ 10 .^ (3:6)
+    reltols = 1 ./ 10 .^ (3:6)
+
+    setups1 = [
+        Dict{Symbol, Any}(:alg => RK4(), :tags => [:rk]),
+    ]
+    setups2 = [
+        Dict{Symbol, Any}(:alg => DP5(), :tags => [:rk]),
+        Dict{Symbol, Any}(:alg => BS3(), :tags => [:rk]),
+    ]
+
+    wp1 = WorkPrecisionSet(prob, abstols, reltols, setups1; dt = 1 / 2^4, numruns = 2)
+    wp2 = WorkPrecisionSet(prob, abstols, reltols, setups2; dt = 1 / 2^4, numruns = 2)
+
+    merged = merge_wp_sets(wp1, wp2)
+    @test length(merged) == 3
+    @test merged.names == ["RK4", "DP5", "BS3"]
+end
diff --git a/lib/DiffEqDevTools/test/timeout_tests.jl b/lib/DiffEqDevTools/test/timeout_tests.jl
new file mode 100644
index 00000000000..a4a42a2b005
--- /dev/null
+++ b/lib/DiffEqDevTools/test/timeout_tests.jl
@@ -0,0 +1,28 @@
+using OrdinaryDiffEq, DiffEqDevTools, Test
+using ODEProblemLibrary: prob_ode_linear
+
+prob = prob_ode_linear
+
+@testset "Timeout Tests" begin
+    abstols = 1 ./ 10 .^ (3:6)
+    reltols = 1 ./ 10 .^ (3:6)
+
+    @testset "timeout=nothing preserves behavior" begin
+        wp = WorkPrecision(prob, RK4(), abstols, reltols;
+            name = "RK4", dt = 1 / 2^4, numruns = 2, timeout = nothing)
+        @test !any(isnan, wp.times)
+    end
+
+    @testset "timeout=0.0 causes NaN" begin
+        wp = @test_warn r"Timeout" WorkPrecision(prob, RK4(), abstols, reltols;
+            name = "RK4", dt = 1 / 2^4, numruns = 2, timeout = 0.0)
+        @test all(isnan, wp.times)
+    end
+
+    @testset "WorkPrecisionSet passes timeout" begin
+        setups = [Dict{Symbol, Any}(:alg => RK4())]
+        wp_set = @test_warn r"Timeout" WorkPrecisionSet(prob, abstols, reltols, setups;
+            dt = 1 / 2^4, numruns = 2, timeout = 0.0)
+        @test all(isnan, wp_set[1].times)
+    end
+end