diff --git a/lib/ModelingToolkitBase/test/analysis_points.jl b/lib/ModelingToolkitBase/test/analysis_points.jl
index c4afca1756..90da216cbe 100644
--- a/lib/ModelingToolkitBase/test/analysis_points.jl
+++ b/lib/ModelingToolkitBase/test/analysis_points.jl
@@ -779,6 +779,312 @@ if @isdefined(ModelingToolkit)
         ]
         @test isapprox(fr, reference_fr)
     end
+
+    @testset "isolate_subsystem" begin
+        @testset "basic plant isolation" begin
+            @named P = FirstOrder(k = 1, T = 1)
+            @named C = Blocks.Gain(k = -1)
+
+            eqs = [connect(C.output, :u, P.input), connect(P.output, :y, C.input)]
+            sys = System(eqs, t, systems = [P, C], name = :cl)
+
+            isolated, input_vars, output_vars = isolate_subsystem(sys, :u, :y)
+
+            @test length(ModelingToolkit.get_systems(isolated)) == 1
+            @test nameof(only(ModelingToolkit.get_systems(isolated))) == :P
+            @test isempty(ModelingToolkit.get_eqs(isolated))
+            @test isequal(only(input_vars), P.input.u)
+            @test isequal(only(output_vars), P.output.u)
+        end
+
+        @testset "external components are removed" begin
+            @named P = FirstOrder(k = 1, T = 1)
+            @named C = Blocks.Gain(k = 1)
+            @named add = Blocks.Add(k2 = -1)
+            @named ref = Step()
+
+            eqs = [
+                connect(ref.output, add.input1)
+                connect(P.output, :y, add.input2)
+                connect(add.output, C.input)
+                connect(C.output, :u, P.input)
+            ]
+            sys = System(eqs, t, systems = [P, C, add, ref], name = :cl)
+
+            isolated, input_vars, output_vars = isolate_subsystem(sys, :u, :y)
+
+            @test Set(nameof.(ModelingToolkit.get_systems(isolated))) == Set([:P])
+            @test isempty(ModelingToolkit.get_eqs(isolated))
+            @test isequal(only(input_vars), P.input.u)
+            @test isequal(only(output_vars), P.output.u)
+        end
+
+        @testset "internal connections are preserved" begin
+            @named P1 = FirstOrder(k = 1, T = 1)
+            @named P2 = FirstOrder(k = 1, T = 1)
+            @named C = Blocks.Gain(k = -1)
+
+            eqs = [
+                connect(C.output, :u, P1.input)
+                connect(P1.output, P2.input)
+                connect(P2.output, :y, C.input)
+            ]
+            sys = System(eqs, t, systems = [P1, P2, C], name = :cl)
+
+            isolated, input_vars, output_vars = isolate_subsystem(sys, :u, :y)
+
+            @test Set(nameof.(ModelingToolkit.get_systems(isolated))) == Set([:P1, :P2])
+            @test length(ModelingToolkit.get_eqs(isolated)) == 1
+            @test Symbolics.value(only(ModelingToolkit.get_eqs(isolated)).rhs) isa Connection
+            @test isequal(only(input_vars), P1.input.u)
+            @test isequal(only(output_vars), P2.output.u)
+        end
+
+        @testset "reachability finds intermediate inside components" begin
+            @named P1 = FirstOrder(k = 1, T = 1)
+            @named P_mid = FirstOrder(k = 1, T = 1)
+            @named P2 = FirstOrder(k = 1, T = 1)
+            @named C = Blocks.Gain(k = -1)
+
+            eqs = [
+                connect(C.output, :u, P1.input)
+                connect(P1.output, P_mid.input)
+                connect(P_mid.output, P2.input)
+                connect(P2.output, :y, C.input)
+            ]
+            sys = System(eqs, t, systems = [P1, P_mid, P2, C], name = :cl)
+
+            isolated, input_vars, output_vars = isolate_subsystem(sys, :u, :y)
+
+            @test Set(nameof.(ModelingToolkit.get_systems(isolated))) == Set([:P1, :P_mid, :P2])
+            @test length(ModelingToolkit.get_eqs(isolated)) == 2
+            @test isequal(only(input_vars), P1.input.u)
+            @test isequal(only(output_vars), P2.output.u)
+        end
+
+        @testset "AnalysisPoint object API" begin
+            @named P = FirstOrder(k = 1, T = 1)
+            @named C = Blocks.Gain(k = -1)
+
+            eqs = [connect(C.output, :u, P.input), connect(P.output, :y, C.input)]
+            sys = System(eqs, t, systems = [P, C], name = :cl)
+
+            isolated, input_vars, output_vars = isolate_subsystem(sys, sys.u, sys.y)
+
+            @test Set(nameof.(ModelingToolkit.get_systems(isolated))) == Set([:P])
+            @test isequal(only(input_vars), P.input.u)
+            @test isequal(only(output_vars), P.output.u)
+        end
+
+        @testset "causal variable connectors" begin
+            @named P = FirstOrder(k = 1, T = 1)
+            @named C = Blocks.Gain(k = -1)
+
+            eqs = [
+                connect(C.output.u, :u, P.input.u)
+                connect(P.output.u, :y, C.input.u)
+            ]
+            sys = System(eqs, t, systems = [P, C], name = :cl)
+
+            isolated, input_vars, output_vars = isolate_subsystem(sys, :u, :y)
+
+            @test Set(nameof.(ModelingToolkit.get_systems(isolated))) == Set([:P])
+            @test isempty(ModelingToolkit.get_eqs(isolated))
+            @test isequal(only(input_vars), P.input.u)
+            @test isequal(only(output_vars), P.output.u)
+        end
+
+        @testset "vector of symbol API" begin
+            @named P = FirstOrder(k = 1, T = 1)
+            @named C = Blocks.Gain(k = -1)
+
+            eqs = [connect(C.output, :u, P.input), connect(P.output, :y, C.input)]
+            sys = System(eqs, t, systems = [P, C], name = :cl)
+
+            isolated, input_vars, output_vars = isolate_subsystem(sys, [:u], [:y])
+
+            @test Set(nameof.(ModelingToolkit.get_systems(isolated))) == Set([:P])
+            @test isequal(only(input_vars), P.input.u)
+            @test isequal(only(output_vars), P.output.u)
+        end
+
+        @testset "nested analysis points" begin
+            @named P = FirstOrder(k = 1, T = 1)
+            @named C = Blocks.Gain(k = -1)
+
+            # APs live inside `inner`, not at the root level
+            inner_eqs = [connect(C.output, :u, P.input), connect(P.output, :y, C.input)]
+            @named inner = System(inner_eqs, t, systems = [P, C])
+            @named root = System(Equation[], t, systems = [inner])
+
+            # Access APs through the nested hierarchy using AnalysisPoint objects
+            isolated, input_vars, output_vars = isolate_subsystem(
+                root, root.inner.u, root.inner.y
+            )
+
+            # root is returned; its only direct child is a trimmed inner containing only P
+            @test Set(nameof.(ModelingToolkit.get_systems(isolated))) == Set([:inner])
+            inner_isolated = only(ModelingToolkit.get_systems(isolated))
+            @test Set(nameof.(ModelingToolkit.get_systems(inner_isolated))) == Set([:P])
+            @test isempty(ModelingToolkit.get_eqs(isolated))
+            @test isempty(ModelingToolkit.get_eqs(inner_isolated))
+            @test isequal(only(input_vars), P.input.u)
+            @test isequal(only(output_vars), P.output.u)
+        end
+
+        @testset "nested analysis points - symbol API" begin
+            @named P = FirstOrder(k = 1, T = 1)
+            @named C = Blocks.Gain(k = -1)
+
+            inner_eqs = [connect(C.output, :u, P.input), connect(P.output, :y, C.input)]
+            @named inner = System(inner_eqs, t, systems = [P, C])
+            @named root = System(Equation[], t, systems = [inner])
+
+            # Access APs by their full namespaced symbol
+            isolated, input_vars, output_vars = isolate_subsystem(
+                root, nameof(inner.u), nameof(inner.y)
+            )
+
+            @test Set(nameof.(ModelingToolkit.get_systems(isolated))) == Set([:inner])
+            inner_isolated = only(ModelingToolkit.get_systems(isolated))
+            @test Set(nameof.(ModelingToolkit.get_systems(inner_isolated))) == Set([:P])
+            @test isequal(only(input_vars), P.input.u)
+            @test isequal(only(output_vars), P.output.u)
+        end
+
+        @testset "nested with external components at outer level" begin
+            @named P = FirstOrder(k = 1, T = 1)
+            @named C = Blocks.Gain(k = -1)
+            @named ref = Step()
+
+            # The APs bounding the plant live inside `inner`
+            inner_eqs = [connect(C.output, :u, P.input), connect(P.output, :y, C.input)]
+            @named inner = System(inner_eqs, t, systems = [P, C])
+
+            # `ref` exists at the outer level — it must not bleed into the isolated result
+            outer_eqs = [connect(ref.output, inner.C.input)]
+            @named root = System(outer_eqs, t, systems = [inner, ref])
+
+            isolated, input_vars, output_vars = isolate_subsystem(
+                root, root.inner.u, root.inner.y
+            )
+
+            # root is returned; ref is stripped, inner is trimmed to only P
+            @test Set(nameof.(ModelingToolkit.get_systems(isolated))) == Set([:inner])
+            inner_isolated = only(ModelingToolkit.get_systems(isolated))
+            @test Set(nameof.(ModelingToolkit.get_systems(inner_isolated))) == Set([:P])
+            @test isempty(ModelingToolkit.get_eqs(isolated))
+            @test isempty(ModelingToolkit.get_eqs(inner_isolated))
+            @test isequal(only(input_vars), P.input.u)
+            @test isequal(only(output_vars), P.output.u)
+        end
+
+        @testset "mixed nesting levels" begin
+            @named P = FirstOrder(k = 1, T = 1)
+            @named C = Blocks.Gain(k = -1)
+            @named A = Step()
+
+            # AP :y lives inside `inner`; AP :u lives at root level
+            inner_eqs = [connect(P.output, :y, C.input)]
+            @named inner = System(inner_eqs, t, systems = [P, C])
+
+            # A drives P.input through AP :u at the root level
+            outer_eqs = [connect(A.output, :u, inner.P.input)]
+            @named root = System(outer_eqs, t, systems = [A, inner])
+
+            # :u is at root level, inner.y is nested — different nesting levels
+            isolated, input_vars, output_vars = isolate_subsystem(
+                root, :u, root.inner.y
+            )
+
+            # root is returned; A is stripped, inner is trimmed to only P (C removed)
+            @test Set(nameof.(ModelingToolkit.get_systems(isolated))) == Set([:inner])
+            inner_isolated = only(ModelingToolkit.get_systems(isolated))
+            @test Set(nameof.(ModelingToolkit.get_systems(inner_isolated))) == Set([:P])
+            @test isempty(ModelingToolkit.get_eqs(isolated))
+            @test isempty(ModelingToolkit.get_eqs(inner_isolated))
+            # input_var: from root-level AP :u, connector is inner.P.input (root-namespaced)
+            @test isequal(only(input_vars), inner.P.input.u)
+            # output_var: from inner-level AP :y, connector is P.output (inner-namespaced)
+            @test isequal(only(output_vars), P.output.u)
+        end
+
+        @testset "deep nesting — isolate middle two of four" begin
+            @named A = Blocks.Gain(k = 1)
+            @named B = FirstOrder(k = 1, T = 1)
+            @named C = FirstOrder(k = 1, T = 2)
+            @named D = Blocks.Gain(k = 1)
+
+            # Four components in series inside `inner`; APs bound B and C (the middle two)
+            inner_eqs = [
+                connect(A.output, :ap_in, B.input),
+                connect(B.output, :bc, C.input),
+                connect(C.output, :ap_out, D.input),
+            ]
+            @named inner = System(inner_eqs, t, systems = [A, B, C, D])
+            @named root = System(Equation[], t, systems = [inner])
+
+            isolated, input_vars, output_vars = isolate_subsystem(
+                root, root.inner.ap_in, root.inner.ap_out
+            )
+
+            # root is returned; inner is trimmed to contain only B and C
+            @test Set(nameof.(ModelingToolkit.get_systems(isolated))) == Set([:inner])
+            inner_isolated = only(ModelingToolkit.get_systems(isolated))
+            @test Set(nameof.(ModelingToolkit.get_systems(inner_isolated))) == Set([:B, :C])
+            # The B→C connection equation is preserved; A and D boundary APs are removed
+            @test length(ModelingToolkit.get_eqs(inner_isolated)) == 1
+            @test isequal(only(input_vars), B.input.u)
+            @test isequal(only(output_vars), C.output.u)
+            # Container levels have no own variables/parameters/observed/defaults
+            @test isempty(ModelingToolkit.get_unknowns(isolated))
+            @test isempty(ModelingToolkit.get_ps(isolated))
+            @test isempty(ModelingToolkit.get_unknowns(inner_isolated))
+            @test isempty(ModelingToolkit.get_ps(inner_isolated))
+        end
+
+        @testset "container-level variables, parameters, and equations are stripped" begin
+            @named P = FirstOrder(k = 1, T = 1)
+            @named Q = FirstOrder(k = 1, T = 2)
+            @named R = FirstOrder(k = 1, T = 3)
+            @named S = Blocks.Gain(k = 1)
+
+            # Declare an extra variable and parameter at the inner (container) level
+            @variables extra_state(t) = 0.0
+            @parameters extra_gain = 2.0
+
+            inner_eqs = [
+                connect(P.output, :ap_in, Q.input),
+                connect(Q.output, :qr, R.input),
+                connect(R.output, :ap_out, S.input),
+                # Plain algebraic equation declared at the container level — must be removed
+                extra_state ~ extra_gain * Q.output.u,
+            ]
+            # inner has its own unknowns, ps, defaults, and a non-connection equation
+            inner = System(
+                inner_eqs, t, [extra_state], [extra_gain];
+                name = :inner, systems = [P, Q, R, S]
+            )
+            @named root = System(Equation[], t, systems = [inner])
+
+            isolated, input_vars, output_vars = isolate_subsystem(
+                root, root.inner.ap_in, root.inner.ap_out
+            )
+
+            inner_isolated = only(ModelingToolkit.get_systems(isolated))
+            @test Set(nameof.(ModelingToolkit.get_systems(inner_isolated))) == Set([:Q, :R])
+            # The Q→R connection AP is kept; extra_state equation and boundary APs are removed
+            @test length(ModelingToolkit.get_eqs(inner_isolated)) == 1
+            # extra_state, extra_gain, their defaults, and the algebraic equation are stripped
+            @test isempty(ModelingToolkit.get_unknowns(inner_isolated))
+            @test isempty(ModelingToolkit.get_ps(inner_isolated))
+            @test isempty(ModelingToolkit.get_observed(inner_isolated))
+            @test isempty(ModelingToolkit.get_initial_conditions(inner_isolated))
+            # root is also clean
+            @test isempty(ModelingToolkit.get_unknowns(isolated))
+            @test isempty(ModelingToolkit.get_ps(isolated))
+        end
+    end
 end
 
 using DynamicQuantities
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index a2f545bce3..71ef030e21 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -163,6 +163,7 @@ export TearingState
 export Clock, SolverStepClock, TimeDomain
 export get_sensitivity_function, get_comp_sensitivity_function,
     get_looptransfer_function, get_sensitivity, get_comp_sensitivity, get_looptransfer
+export isolate_subsystem
 
 function FMIComponent end
 
diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl
index 769c56a6e0..6e1399a622 100644
--- a/src/systems/analysis_points.jl
+++ b/src/systems/analysis_points.jl
@@ -182,6 +182,220 @@ function linearization_function(
     return linearization_function(system_modifier(sys), input_vars, output_vars; kwargs...)
 end
 
+"""
+    sys, input_vars, output_vars = isolate_subsystem(sys, input_aps, output_aps)
+
+Isolate the part of the unsimplified, hierarchical system `sys` bounded by the input
+analysis points `input_aps` and the output analysis points `output_aps`. The returned
+`sys` contains only the subsystems between the boundary analysis points at every level
+of the hierarchy; all upstream and downstream components, and all equations involving
+them, are removed.
+
+Boundary analysis points may reside at any level of the hierarchy and in different
+branches of the subsystem tree.
+
+Returns:
+- `sys`: The system with only the isolated subsystems and their internal connections.
+- `input_vars`: Variables at the inside face of each input analysis point.
+- `output_vars`: Variables at the inside face of each output analysis point.
+"""
+function isolate_subsystem(
+        sys::AbstractSystem,
+        input_aps::Union{Symbol, Vector{Symbol}, AnalysisPoint, Vector{AnalysisPoint}},
+        output_aps::Union{Symbol, Vector{Symbol}, AnalysisPoint, Vector{AnalysisPoint}}
+    )
+    input_aps = canonicalize_ap(sys, input_aps)
+    output_aps = canonicalize_ap(sys, output_aps)
+
+    # Step 1: index every named subsystem in the hierarchy.
+    # paths[i] is the path of names from sys to that subsystem, e.g. [:inner, :P].
+    paths = Vector{Symbol}[]
+    path_to_idx = Dict{Vector{Symbol}, Int}()
+    function _collect_subsystems!(cur, path)
+        for s in get_systems(cur)
+            p = [path; nameof(s)]
+            push!(paths, p)
+            path_to_idx[p] = length(paths)
+            _collect_subsystems!(s, p)
+        end
+        return
+    end
+    _collect_subsystems!(sys, Symbol[])
+
+    g = SimpleGraph(length(paths))
+
+    input_ap_names = Set{Symbol}(nameof(ap) for ap in input_aps)
+    output_ap_names = Set{Symbol}(nameof(ap) for ap in output_aps)
+    boundary_ap_names = union(input_ap_names, output_ap_names)
+    sys_root_name = nameof(sys)
+
+    # Reconstruct the full canonical AP name (as produced by canonicalize_ap) from the
+    # local name seen in an equation at depth `parent_path` within `sys`.
+    function _full_ap_name(parent_path, local_name)
+        return Symbol(join(string.([sys_root_name; parent_path; local_name]), NAMESPACE_SEPARATOR))
+    end
+
+    # Map a connector (port System or signal SymbolicT) at a given parent path to the
+    # component path by stripping the last namespace segment (port or variable name).
+    function _conn_to_path(conn, parent_path)
+        segs = conn isa AbstractSystem ? namespace_hierarchy(nameof(conn)) :
+            namespace_hierarchy(getname(unwrap(conn)))
+        length(segs) <= 1 && return nothing
+        return [parent_path; segs[1:(end - 1)]]
+    end
+
+    function _try_add_edge!(p1, p2)
+        i1 = get(path_to_idx, p1, nothing)
+        i2 = get(path_to_idx, p2, nothing)
+        (i1 === nothing || i2 === nothing || i1 == i2) && return
+        return add_edge!(g, i1, i2)
+    end
+
+    # Step 2: walk every equation at every level of the hierarchy.
+    # For boundary AP equations, cut the edge and record which side is inside.
+    # For all other equations, add edges between the connected components.
+    inside_seeds = Int[]
+    input_vars = SymbolicT[]
+    output_vars = SymbolicT[]
+
+    function _build_graph!(cur, parent_path)
+        for eq in get_eqs(cur)
+            lhs_val = value(eq.lhs)
+            rhs_val = value(eq.rhs)
+            if lhs_val isa AnalysisPoint
+                ap_data = rhs_val::AnalysisPoint
+                fname = _full_ap_name(parent_path, nameof(ap_data))
+                in_conn = ap_data.input
+                out_conns = something(ap_data.outputs, [])
+                if fname in boundary_ap_names
+                    if fname in input_ap_names
+                        for c in out_conns
+                            push!(input_vars, ap_var(c))
+                            p = _conn_to_path(c, parent_path)
+                            if p !== nothing
+                                idx = get(path_to_idx, p, nothing)
+                                idx !== nothing && push!(inside_seeds, idx)
+                            end
+                        end
+                    end
+                    if fname in output_ap_names
+                        if in_conn !== nothing
+                            push!(output_vars, ap_var(in_conn))
+                            p = _conn_to_path(in_conn, parent_path)
+                            if p !== nothing
+                                idx = get(path_to_idx, p, nothing)
+                                idx !== nothing && push!(inside_seeds, idx)
+                            end
+                        end
+                    end
+                    # Do not add a graph edge — this connection is cut.
+                else
+                    in_path = in_conn !== nothing ? _conn_to_path(in_conn, parent_path) : nothing
+                    for c in out_conns
+                        out_path = _conn_to_path(c, parent_path)
+                        out_path === nothing && continue
+                        in_path !== nothing && _try_add_edge!(in_path, out_path)
+                    end
+                end
+            elseif rhs_val isa Connection
+                conn_list = get_systems(rhs_val)
+                conn_list === nothing && continue
+                cps = [
+                    p for c in conn_list
+                        for p in (_conn_to_path(c, parent_path),) if p !== nothing
+                ]
+                for i in eachindex(cps), j in (i + 1):length(cps)
+                    _try_add_edge!(cps[i], cps[j])
+                end
+            end
+        end
+        for s in get_systems(cur)
+            _build_graph!(s, [parent_path; nameof(s)])
+        end
+        return
+    end
+    _build_graph!(sys, Symbol[])
+
+    # Step 3: BFS from inside seeds to find all reachable (inside) components.
+    inside = falses(length(paths))
+    queue = unique!(copy(inside_seeds))
+    for v in queue
+        inside[v] = true
+    end
+    qi = 1
+    while qi <= length(queue)
+        v = queue[qi]
+        qi += 1
+        for w in Graphs.neighbors(g, v)
+            if !inside[w]
+                inside[w] = true
+                push!(queue, w)
+            end
+        end
+    end
+
+    # Build a prefix set so that has_inside(path) is an O(1) lookup:
+    # a path is "has inside" if any inside component lives under it.
+    inside_prefixes = Set{Vector{Symbol}}()
+    for i in eachindex(paths)
+        inside[i] || continue
+        p = paths[i]
+        for k in eachindex(p)
+            push!(inside_prefixes, p[1:k])
+        end
+    end
+    has_inside(path) = path in inside_prefixes
+
+    # Step 4: rebuild the system hierarchy keeping only inside components and equations.
+    # Systems that are directly inside (their path is in the `inside` set) are returned
+    # as-is — their internal structure belongs to the isolated subsystem.
+    # Systems that are containers (in inside_prefixes only because they wrap inside
+    # components) have their equations filtered and their own variables/parameters cleared.
+    function _reconstruct!(cur, parent_path)
+        idx = get(path_to_idx, parent_path, nothing)
+        if idx !== nothing && inside[idx]
+            # Directly-inside component — preserve it entirely.
+            return cur
+        end
+
+        # Container: filter subsystems and equations, clear own vars/params so that
+        # nothing from outside the isolated region leaks into the result.
+        new_systems = [
+            _reconstruct!(s, [parent_path; nameof(s)])
+                for s in get_systems(cur) if has_inside([parent_path; nameof(s)])
+        ]
+        new_eqs = filter(get_eqs(cur)) do eq
+            lhs_val = value(eq.lhs)
+            rhs_val = value(eq.rhs)
+            if lhs_val isa AnalysisPoint
+                ap_data = rhs_val::AnalysisPoint
+                _full_ap_name(parent_path, nameof(ap_data)) in boundary_ap_names && return false
+                in_conn = ap_data.input
+                out_conns = something(ap_data.outputs, [])
+                all_conns = in_conn === nothing ? out_conns : [in_conn; out_conns]
+                all_paths = [
+                    p for c in all_conns
+                        for p in (_conn_to_path(c, parent_path),) if p !== nothing
+                ]
+                return all(has_inside, all_paths)
+            elseif rhs_val isa Connection
+                conn_list = get_systems(rhs_val)
+                conn_list === nothing && return true
+                cps = [
+                    p for c in conn_list
+                        for p in (_conn_to_path(c, parent_path),) if p !== nothing
+                ]
+                return all(has_inside, cps)
+            else
+                return false
+            end
+        end
+        return System(new_eqs, get_iv(cur), SymbolicT[], SymbolicT[]; name = nameof(cur), systems = new_systems)
+    end
+
+    return _reconstruct!(sys, Symbol[]), input_vars, output_vars
+end
+
 @doc """
     get_sensitivity(sys, ap::AnalysisPoint; kwargs)
     get_sensitivity(sys, ap_name::Symbol; kwargs)