diff --git a/Project.toml b/Project.toml
index 3c6126b7c0..859eb27170 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "9.59.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -73,11 +74,12 @@ MTKBifurcationKitExt = "BifurcationKit"
 MTKChainRulesCoreExt = "ChainRulesCore"
 MTKDeepDiffsExt = "DeepDiffs"
 MTKHomotopyContinuationExt = "HomotopyContinuation"
-MTKLabelledArraysExt = "LabelledArrays"
 MTKInfiniteOptExt = "InfiniteOpt"
+MTKLabelledArraysExt = "LabelledArrays"
 
 [compat]
 AbstractTrees = "0.3, 0.4"
+Accessors = "0.1.36"
 ArrayInterface = "6, 7"
 BifurcationKit = "0.4"
 BlockArrays = "1.1"
diff --git a/docs/Project.toml b/docs/Project.toml
index a358455503..29c07f5c55 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
@@ -23,6 +24,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
+Accessors = "0.1.36"
 BenchmarkTools = "1.3"
 BifurcationKit = "0.4"
 DataInterpolations = "6.5"
diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md
index 23e1e6d7d1..6d3a90a8d4 100644
--- a/docs/src/basics/Events.md
+++ b/docs/src/basics/Events.md
@@ -412,16 +412,16 @@ is below `furnace_on_threshold` and off when above `furnace_off_threshold`, whil
 in between. To do this, we create two continuous callbacks:
 
 ```@example events
-using Setfield
+using Accessors
 furnace_disable = ModelingToolkit.SymbolicContinuousCallback(
     [temp ~ furnace_off_threshold],
     ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i
-        @set! x.furnace_on = false
+        @reset x.furnace_on = false
     end)
 furnace_enable = ModelingToolkit.SymbolicContinuousCallback(
     [temp ~ furnace_on_threshold],
     ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i
-        @set! x.furnace_on = true
+        @reset x.furnace_on = true
     end)
 ```
 
@@ -432,7 +432,7 @@ You can also write
 ```julia
 [temp ~ furnace_off_threshold] => ModelingToolkit.ImperativeAffect(modified = (;
     furnace_on)) do x, o, i, c
-    @set! x.furnace_on = false
+    @reset x.furnace_on = false
 end
 ```
 
@@ -462,7 +462,7 @@ f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple
 The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions.
 In our example, if `furnace_on` is `false`, then the value of the `x` that's passed in as `modified` will be `(furnace_on = false)`.
 The modified values should be passed out in the same format: to set `furnace_on` to `true` we need to return a tuple `(furnace_on = true)`.
-The examples does this with Setfield, recreating the result tuple before returning it; the returned tuple may optionally be missing values as
+The examples does this with Accessors, recreating the result tuple before returning it; the returned tuple may optionally be missing values as
 well, in which case those values will not be written back to the problem.
 
 Accordingly, we can now interpret the `ImperativeAffect` definitions to mean that when `temp = furnace_off_threshold` we
@@ -542,18 +542,18 @@ In our encoder, we interpret this as occlusion or nonocclusion of the sensor, up
 ```@example events
 qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0],
     ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, c, i
-        @set! x.hA = x.qA
-        @set! x.hB = o.qB
-        @set! x.qA = 1
-        @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
+        @reset x.hA = x.qA
+        @reset x.hB = o.qB
+        @reset x.qA = 1
+        @reset x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
         x
     end,
     affect_neg = ModelingToolkit.ImperativeAffect(
         (; qA, hA, hB, cnt), (; qB)) do x, o, c, i
-        @set! x.hA = x.qA
-        @set! x.hB = o.qB
-        @set! x.qA = 0
-        @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
+        @reset x.hA = x.qA
+        @reset x.hB = o.qB
+        @reset x.qA = 0
+        @reset x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
         x
     end)
 ```
@@ -566,10 +566,10 @@ Instead, we can use right root finding:
 ```@example events
 qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0],
     ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA, theta)) do x, o, c, i
-        @set! x.hA = o.qA
-        @set! x.hB = x.qB
-        @set! x.qB = clamp(sign(cos(100 * o.theta - π / 2)), 0.0, 1.0)
-        @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB)
+        @reset x.hA = o.qA
+        @reset x.hB = x.qB
+        @reset x.qB = clamp(sign(cos(100 * o.theta - π / 2)), 0.0, 1.0)
+        @reset x.cnt += decoder(x.hA, x.hB, o.qA, x.qB)
         x
     end; rootfind = SciMLBase.RightRootFind)
 ```
diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
index e5b57efc5a..2f42815630 100644
--- a/ext/MTKBifurcationKitExt.jl
+++ b/ext/MTKBifurcationKitExt.jl
@@ -3,7 +3,7 @@ module MTKBifurcationKitExt
 ### Preparations ###
 
 # Imports
-using ModelingToolkit, Setfield
+using ModelingToolkit, Accessors
 import BifurcationKit
 
 ### Observable Plotting Handling ###
@@ -94,7 +94,7 @@ function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem,
     if !ModelingToolkit.iscomplete(nsys)
         error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `BifurcationProblem`")
     end
-    @set! nsys.index_cache = nothing # force usage of a parameter vector instead of `MTKParameters`
+    @reset nsys.index_cache = nothing # force usage of a parameter vector instead of `MTKParameters`
     # Creates F and J functions.
     ofun = NonlinearFunction(nsys; jac = jac)
     F = ofun.f
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 10aba4d8a9..dea49dee92 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -31,7 +31,8 @@ using JumpProcesses
 using DataStructures
 using Base.Threads
 using Latexify, Unitful, ArrayInterface
-using Setfield, ConstructionBase
+import Setfield
+using Accessors, ConstructionBase
 import Libdl
 using DocStringExtensions
 using Base: RefValue
diff --git a/src/bipartite_graph.jl b/src/bipartite_graph.jl
index b6665646c9..8db3c7d803 100644
--- a/src/bipartite_graph.jl
+++ b/src/bipartite_graph.jl
@@ -14,7 +14,7 @@ using DocStringExtensions
 using UnPack
 using SparseArrays
 using Graphs
-using Setfield
+using Accessors
 
 ### Matching
 struct Unassigned
diff --git a/src/inputoutput.jl b/src/inputoutput.jl
index 6bdcac6dd4..3ceb689cfb 100644
--- a/src/inputoutput.jl
+++ b/src/inputoutput.jl
@@ -315,12 +315,12 @@ function inputs_to_parameters!(state::TransformationState, io)
         @assert new_v > 0
         new_var_to_diff[new_i] = new_v
     end
-    @set! structure.var_to_diff = complete(new_var_to_diff)
-    @set! structure.graph = complete(new_graph)
+    @reset structure.var_to_diff = complete(new_var_to_diff)
+    @reset structure.graph = complete(new_graph)
 
-    @set! sys.eqs = isempty(input_to_parameters) ? equations(sys) :
+    @reset sys.eqs = isempty(input_to_parameters) ? equations(sys) :
                     fast_substitute(equations(sys), input_to_parameters)
-    @set! sys.unknowns = setdiff(unknowns(sys), keys(input_to_parameters))
+    @reset sys.unknowns = setdiff(unknowns(sys), keys(input_to_parameters))
     ps = parameters(sys)
 
     if io !== nothing
@@ -334,11 +334,11 @@ function inputs_to_parameters!(state::TransformationState, io)
         new_parameters = new_parameters[permutation]
     end
 
-    @set! sys.ps = [ps; new_parameters]
+    @reset sys.ps = [ps; new_parameters]
 
-    @set! state.sys = sys
-    @set! state.fullvars = new_fullvars
-    @set! state.structure = structure
+    @reset state.sys = sys
+    @reset state.fullvars = new_fullvars
+    @reset state.structure = structure
     base_params = length(ps)
     return state, (base_params + 1):(base_params + length(new_parameters)) # (1:length(new_parameters)) .+ base_params
 end
diff --git a/src/structural_transformation/StructuralTransformations.jl b/src/structural_transformation/StructuralTransformations.jl
index 1220d517cc..50189c04dc 100644
--- a/src/structural_transformation/StructuralTransformations.jl
+++ b/src/structural_transformation/StructuralTransformations.jl
@@ -1,6 +1,6 @@
 module StructuralTransformations
 
-using Setfield: @set!, @set
+using Accessors: @set, @reset
 using UnPack: @unpack
 
 using Symbolics: unwrap, linear_expansion, fast_substitute
diff --git a/src/structural_transformation/codegen.jl b/src/structural_transformation/codegen.jl
index d2ca5b8748..97cdd47f0d 100644
--- a/src/structural_transformation/codegen.jl
+++ b/src/structural_transformation/codegen.jl
@@ -300,7 +300,7 @@ function build_torn_function(sys;
         rhss)
 
     unknown_vars = Any[fullvars[i] for i in unknowns_idxs]
-    @set! sys.solved_unknowns = unknown_vars
+    @reset sys.solved_unknowns = unknown_vars
 
     pre = get_postprocess_fbody(sys)
     cpre = get_preprocess_constants(rhss)
diff --git a/src/structural_transformation/pantelides.jl b/src/structural_transformation/pantelides.jl
index b6877d65f8..1567186cc1 100644
--- a/src/structural_transformation/pantelides.jl
+++ b/src/structural_transformation/pantelides.jl
@@ -65,8 +65,8 @@ function pantelides_reassemble(state::TearingState, var_eq_matching)
         filter(x -> value(x.lhs) !== nothing,
             out_eqs[sort(filter(x -> x !== unassigned, var_eq_matching))]))
 
-    @set! sys.eqs = final_eqs
-    @set! sys.unknowns = final_vars
+    @reset sys.eqs = final_eqs
+    @reset sys.unknowns = final_vars
     return sys
 end
 
diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
index 8161a9572b..eaf2e045b3 100644
--- a/src/structural_transformation/symbolics_tearing.jl
+++ b/src/structural_transformation/symbolics_tearing.jl
@@ -133,9 +133,9 @@ end
 
 function tearing_substitution(sys::AbstractSystem; kwargs...)
     neweqs = full_equations(sys::AbstractSystem; kwargs...)
-    @set! sys.eqs = neweqs
-    @set! sys.substitutions = nothing
-    @set! sys.schedule = nothing
+    @reset sys.eqs = neweqs
+    @reset sys.substitutions = nothing
+    @reset sys.schedule = nothing
 end
 
 function tearing_assignments(sys::AbstractSystem)
@@ -563,10 +563,10 @@ function tearing_reassemble(state::TearingState, var_eq_matching,
     diff_to_var = invview(var_to_diff)
 
     old_fullvars = fullvars
-    @set! state.structure.graph = complete(graph)
-    @set! state.structure.var_to_diff = var_to_diff
-    @set! state.structure.eq_to_diff = eq_to_diff
-    @set! state.fullvars = fullvars = fullvars[invvarsperm]
+    @reset state.structure.graph = complete(graph)
+    @reset state.structure.var_to_diff = var_to_diff
+    @reset state.structure.eq_to_diff = eq_to_diff
+    @reset state.fullvars = fullvars = fullvars[invvarsperm]
     ispresent = let var_to_diff = var_to_diff, graph = graph
         i -> (!isempty(𝑑neighbors(graph, i)) ||
               (var_to_diff[i] !== nothing && !isempty(𝑑neighbors(graph, var_to_diff[i]))))
@@ -590,24 +590,24 @@ function tearing_reassemble(state::TearingState, var_eq_matching,
             push!(unknowns, old_fullvars[v])
         end
     end
-    @set! sys.unknowns = unknowns
+    @reset sys.unknowns = unknowns
 
     obs, subeqs, deps = cse_and_array_hacks(
         obs, subeqs, unknowns, neweqs; cse = cse_hack, array = array_hack)
 
-    @set! sys.eqs = neweqs
-    @set! sys.observed = obs
+    @reset sys.eqs = neweqs
+    @reset sys.observed = obs
 
-    @set! sys.substitutions = Substitutions(subeqs, deps)
+    @reset sys.substitutions = Substitutions(subeqs, deps)
 
     # Only makes sense for time-dependent
     # TODO: generalize to SDE
     if sys isa ODESystem
-        @set! sys.schedule = Schedule(var_eq_matching, dummy_sub)
+        @reset sys.schedule = Schedule(var_eq_matching, dummy_sub)
     end
     sys = schedule(sys)
-    @set! state.sys = sys
-    @set! sys.tearing_state = state
+    @reset state.sys = sys
+    @reset sys.tearing_state = state
     return invalidate_cache!(sys)
 end
 
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 168260ae69..0c6368390d 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -859,7 +859,7 @@ is scheduled after tearing based simplifications where equations are converted
 into assignments.
 """
 function schedule(sys::AbstractSystem)
-    has_schedule(sys) ? sys : (@set! sys.isscheduled = true)
+    has_schedule(sys) ? sys : (@reset sys.isscheduled = true)
 end
 
 """
@@ -896,7 +896,7 @@ function complete(sys::AbstractSystem; split = true, flatten = true)
     collect_scoped_vars!(newunknowns, newparams, sys, iv; depth = -1)
     # don't update unknowns to not disturb `structural_simplify` order
     # `GlobalScope`d unknowns will be picked up and added there
-    @set! sys.ps = unique!(vcat(get_ps(sys), collect(newparams)))
+    @reset sys.ps = unique!(vcat(get_ps(sys), collect(newparams)))
 
     if flatten
         eqs = equations(sys)
@@ -907,12 +907,12 @@ function complete(sys::AbstractSystem; split = true, flatten = true)
         end
         newsys = ModelingToolkit.flatten(newsys)
         if has_parent(newsys) && get_parent(sys) === nothing
-            @set! newsys.parent = complete(sys; split = false, flatten = false)
+            @reset newsys.parent = complete(sys; split = false, flatten = false)
         end
         sys = newsys
     end
     if split && has_index_cache(sys)
-        @set! sys.index_cache = IndexCache(sys)
+        @reset sys.index_cache = IndexCache(sys)
         all_ps = get_ps(sys)
         if !isempty(all_ps)
             # reorder parameters by portions
@@ -949,15 +949,15 @@ function complete(sys::AbstractSystem; split = true, flatten = true)
                 ordered_ps = vcat(
                     ordered_ps, reduce(vcat, ps_split[2:end]; init = eltype(ordered_ps)[]))
             end
-            @set! sys.ps = ordered_ps
+            @reset sys.ps = ordered_ps
         end
     elseif has_index_cache(sys)
-        @set! sys.index_cache = nothing
+        @reset sys.index_cache = nothing
     end
     if isdefined(sys, :initializesystem) && get_initializesystem(sys) !== nothing
-        @set! sys.initializesystem = complete(get_initializesystem(sys); split)
+        @reset sys.initializesystem = complete(get_initializesystem(sys); split)
     end
-    isdefined(sys, :complete) ? (@set! sys.complete = true) : sys
+    isdefined(sys, :complete) ? (@reset sys.complete = true) : sys
 end
 
 for prop in [:eqs
@@ -1062,6 +1062,7 @@ end
 function Setfield.get(obj::AbstractSystem, ::Setfield.PropertyLens{field}) where {field}
     getfield(obj, field)
 end
+
 @generated function ConstructionBase.setproperties(obj::AbstractSystem, patch::NamedTuple)
     if issubset(fieldnames(patch), fieldnames(obj))
         args = map(fieldnames(obj)) do fn
@@ -2229,9 +2230,9 @@ function component_post_processing(expr, isconnector)
             if $isdefined(res, :gui_metadata) && $getfield(res, :gui_metadata) === nothing
                 name = $(Meta.quot(fname))
                 if $isconnector
-                    $Setfield.@set!(res.connector_type=$connector_type(res))
+                    $Accessors.@reset(res.connector_type=$connector_type(res))
                 end
-                $Setfield.@set!(res.gui_metadata=$GUIMetadata($GlobalRef(@__MODULE__, name)))
+                $Accessors.@reset(res.gui_metadata=$GUIMetadata($GlobalRef(@__MODULE__, name)))
             else
                 res
             end
@@ -2287,10 +2288,10 @@ function debug_system(sys::AbstractSystem)
         error("debug_system only works on systems with no sub-systems!")
     end
     if has_eqs(sys)
-        @set! sys.eqs = debug_sub.(equations(sys))
+        @reset sys.eqs = debug_sub.(equations(sys))
     end
     if has_observed(sys)
-        @set! sys.observed = debug_sub.(observed(sys))
+        @reset sys.observed = debug_sub.(observed(sys))
     end
     return sys
 end
@@ -2301,7 +2302,7 @@ function eliminate_constants(sys::AbstractSystem)
         eq_cs = collect_constants(eqs)
         if !isempty(eq_cs)
             new_eqs = eliminate_constants(eqs, eq_cs)
-            @set! sys.eqs = new_eqs
+            @reset sys.eqs = new_eqs
         end
     end
     return sys
@@ -2372,7 +2373,7 @@ function linearization_function(sys::AbstractSystem, inputs,
     if zero_dummy_der
         dummyder = setdiff(unknowns(ssys), unknowns(sys))
         defs = Dict(x => 0.0 for x in dummyder)
-        @set! ssys.defaults = merge(defs, defaults(ssys))
+        @reset ssys.defaults = merge(defs, defaults(ssys))
         op = merge(defs, op)
     end
     sys = ssys
@@ -2392,7 +2393,7 @@ function linearization_function(sys::AbstractSystem, inputs,
             @warn "Initialization system is underdetermined. No equations for $(missing_unknowns). Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false."
         end
         new_parameters = [parameters(initsys); missing_unknowns]
-        @set! initsys.ps = new_parameters
+        @reset initsys.ps = new_parameters
         initsys = complete(initsys)
     end
 
@@ -3051,10 +3052,10 @@ See also [`extend`](@ref).
 function compose(sys::AbstractSystem, systems::AbstractArray; name = nameof(sys))
     nsys = length(systems)
     nsys == 0 && return sys
-    @set! sys.name = name
-    @set! sys.systems = [get_systems(sys); systems]
+    @reset sys.name = name
+    @reset sys.systems = [get_systems(sys); systems]
     if has_is_dde(sys)
-        @set! sys.is_dde = _check_if_dde(equations(sys), get_iv(sys), get_systems(sys))
+        @reset sys.is_dde = _check_if_dde(equations(sys), get_iv(sys), get_systems(sys))
     end
     newunknowns = OrderedSet()
     newparams = OrderedSet()
@@ -3062,8 +3063,8 @@ function compose(sys::AbstractSystem, systems::AbstractArray; name = nameof(sys)
     for ssys in systems
         collect_scoped_vars!(newunknowns, newparams, ssys, iv)
     end
-    @set! sys.unknowns = unique!(vcat(get_unknowns(sys), collect(newunknowns)))
-    @set! sys.ps = unique!(vcat(get_ps(sys), collect(newparams)))
+    @reset sys.unknowns = unique!(vcat(get_unknowns(sys), collect(newunknowns)))
+    @reset sys.ps = unique!(vcat(get_ps(sys), collect(newparams)))
     return sys
 end
 function compose(syss...; name = nameof(first(syss)))
@@ -3116,7 +3117,7 @@ function Symbolics.substitute(sys::AbstractSystem, rules::Union{Vector{<:Pair},
         dict = todict(rules)
         systems = get_systems(sys)
         # post-walk to avoid infinite recursion
-        @set! sys.systems = map(Base.Fix2(substitute, dict), systems)
+        @reset sys.systems = map(Base.Fix2(substitute, dict), systems)
         something(get(rules, nameof(sys), nothing), sys)
     elseif sys isa ODESystem
         rules = todict(map(r -> Symbolics.unwrap(r[1]) => Symbolics.unwrap(r[2]),
diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index fb4fedc920..de3a0e3190 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -99,14 +99,14 @@ function alias_elimination!(state::TearingState; kwargs...)
         eq = eqs[ieq]
         eqs[ieq] = fast_substitute(eq, subs)
     end
-    @set! mm.nparentrows = nsrcs(graph)
-    @set! mm.row_cols = eltype(mm.row_cols)[mm.row_cols[i]
+    @reset mm.nparentrows = nsrcs(graph)
+    @reset mm.row_cols = eltype(mm.row_cols)[mm.row_cols[i]
                                             for (i, eq) in enumerate(mm.nzrows)
                                             if old_to_new_eq[eq] > 0]
-    @set! mm.row_vals = eltype(mm.row_vals)[mm.row_vals[i]
+    @reset mm.row_vals = eltype(mm.row_vals)[mm.row_vals[i]
                                             for (i, eq) in enumerate(mm.nzrows)
                                             if old_to_new_eq[eq] > 0]
-    @set! mm.nzrows = Int[old_to_new_eq[eq] for eq in mm.nzrows if old_to_new_eq[eq] > 0]
+    @reset mm.nzrows = Int[old_to_new_eq[eq] for eq in mm.nzrows if old_to_new_eq[eq] > 0]
 
     for old_ieq in to_expand
         ieq = old_to_new_eq[old_ieq]
@@ -136,7 +136,7 @@ function alias_elimination!(state::TearingState; kwargs...)
     state.structure.var_to_diff = new_var_to_diff
 
     sys = state.sys
-    @set! sys.eqs = eqs
+    @reset sys.eqs = eqs
     state.sys = sys
     return invalidate_cache!(sys), mm
 end
@@ -359,7 +359,7 @@ function alias_eliminate_graph!(state::TransformationState, ils::SparseMatrixCLI
     rk1vars = BitSet(@view pivots[1:rank1])
     for v in solvable_variables
         v in rk1vars && continue
-        @set! ils.nparentrows += 1
+        @reset ils.nparentrows += 1
         push!(ils.nzrows, ils.nparentrows)
         push!(ils.row_cols, [v])
         push!(ils.row_vals, [convert(eltype(ils), 1)])
diff --git a/src/systems/clock_inference.jl b/src/systems/clock_inference.jl
index a92b2aa67c..b2d14e65a7 100644
--- a/src/systems/clock_inference.jl
+++ b/src/systems/clock_inference.jl
@@ -39,8 +39,8 @@ function substitute_sample_time(ci::ClockInference, ts::TearingState)
         end
         eqs[i] = neweq
     end
-    @set! ts.sys.eqs = eqs
-    @set! ci.ts = ts
+    @reset ts.sys.eqs = eqs
+    @reset ci.ts = ts
 end
 
 function substitute_sample_time(eq::Equation, dt)
@@ -182,7 +182,7 @@ function split_system(ci::ClockInference{S}) where {S}
         ts_i = system_subset(ts, ieqs)
         if id != continuous_id
             ts_i = shift_discrete_system(ts_i)
-            @set! ts_i.structure.only_discrete = true
+            @reset ts_i.structure.only_discrete = true
         end
         tss[id] = ts_i
     end
diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl
index 227b4624bf..c6b732e35d 100644
--- a/src/systems/connectors.jl
+++ b/src/systems/connectors.jl
@@ -370,14 +370,14 @@ function generate_connection_set!(connectionsets, domain_csets,
 
     # pre order traversal
     if !isempty(extra_unknowns)
-        @set! sys.unknowns = [get_unknowns(sys); extra_unknowns]
+        @reset sys.unknowns = [get_unknowns(sys); extra_unknowns]
     end
-    @set! sys.systems = map(
+    @reset sys.systems = map(
         s -> generate_connection_set!(connectionsets, domain_csets, s,
             find, replace, scalarize,
             renamespace(namespace, s)),
         subsys)
-    @set! sys.eqs = eqs
+    @reset sys.eqs = eqs
 end
 
 function Base.merge(csets::AbstractVector{<:ConnectionSet}, allouter = false)
@@ -485,9 +485,9 @@ function expand_connections(sys::AbstractSystem, find = nothing, replace = nothi
     ceqs, instream_csets = generate_connection_equations_and_stream_connections(csets)
     _sys = expand_instream(instream_csets, sys; debug = debug, tol = tol)
     sys = flatten(sys, true)
-    @set! sys.eqs = [equations(_sys); ceqs]
+    @reset sys.eqs = [equations(_sys); ceqs]
     d_defs = domain_defaults(sys, domain_csets)
-    @set! sys.defaults = merge(get_defaults(sys), d_defs)
+    @reset sys.defaults = merge(get_defaults(sys), d_defs)
 end
 
 function unnamespace(root, namespace)
@@ -508,7 +508,7 @@ function expand_instream(csets::AbstractVector{<:ConnectionSet}, sys::AbstractSy
         tol = 1e-8)
     subsys = get_systems(sys)
     # post order traversal
-    @set! sys.systems = map(
+    @reset sys.systems = map(
         s -> expand_instream(csets, s,
             renamespace(namespace, nameof(s)),
             namespace; debug, tol),
@@ -672,8 +672,8 @@ function expand_instream(csets::AbstractVector{<:ConnectionSet}, sys::AbstractSy
         println("======================================")
     end
 
-    @set! sys.systems = []
-    @set! sys.eqs = [get_eqs(sys); eqs; subed_eqs; additional_eqs]
+    @reset sys.systems = []
+    @reset sys.eqs = [get_eqs(sys); eqs; subed_eqs; additional_eqs]
     sys
 end
 
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index f4e29346ff..07dc0c663f 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -452,7 +452,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem,
         nothing
     end
 
-    @set! sys.split_idxs = split_idxs
+    @reset sys.split_idxs = split_idxs
 
     ODEFunction{iip, specialize}(f;
         sys = sys,
diff --git a/src/systems/diffeqs/first_order_transform.jl b/src/systems/diffeqs/first_order_transform.jl
index 97fd6460d9..1ce479fbd7 100644
--- a/src/systems/diffeqs/first_order_transform.jl
+++ b/src/systems/diffeqs/first_order_transform.jl
@@ -7,16 +7,16 @@ form by defining new variables which represent the N-1 derivatives.
 function ode_order_lowering(sys::ODESystem)
     iv = get_iv(sys)
     eqs_lowered, new_vars = ode_order_lowering(equations(sys), iv, unknowns(sys))
-    @set! sys.eqs = eqs_lowered
-    @set! sys.unknowns = new_vars
+    @reset sys.eqs = eqs_lowered
+    @reset sys.unknowns = new_vars
     return sys
 end
 
 function dae_order_lowering(sys::ODESystem)
     iv = get_iv(sys)
     eqs_lowered, new_vars = dae_order_lowering(equations(sys), iv, unknowns(sys))
-    @set! sys.eqs = eqs_lowered
-    @set! sys.unknowns = new_vars
+    @reset sys.eqs = eqs_lowered
+    @reset sys.unknowns = new_vars
     return sys
 end
 
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index 34003f40c2..2daca8e2a8 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -766,9 +766,9 @@ function add_accumulations(sys::ODESystem, vars::Vector{<:Pair})
         error("$ints already exist in the system!")
     end
     D = Differential(get_iv(sys))
-    @set! sys.eqs = [eqs; Equation[D(a) ~ v[2] for (a, v) in zip(avars, vars)]]
-    @set! sys.unknowns = [get_unknowns(sys); avars]
-    @set! sys.defaults = merge(get_defaults(sys), Dict(a => 0.0 for a in avars))
+    @reset sys.eqs = [eqs; Equation[D(a) ~ v[2] for (a, v) in zip(avars, vars)]]
+    @reset sys.unknowns = [get_unknowns(sys); avars]
+    @reset sys.defaults = merge(get_defaults(sys), Dict(a => 0.0 for a in avars))
 end
 
 function Base.show(io::IO, mime::MIME"text/plain", sys::ODESystem; hint = true, bold = true)
diff --git a/src/systems/if_lifting.jl b/src/systems/if_lifting.jl
index 53e1ca4957..f20a2977df 100644
--- a/src/systems/if_lifting.jl
+++ b/src/systems/if_lifting.jl
@@ -501,11 +501,11 @@ function IfLifting(sys::ODESystem)
         push!(new_ps, var)
     end
 
-    @set! sys.defaults = new_defaults
-    @set! sys.eqs = eqs
+    @reset sys.defaults = new_defaults
+    @reset sys.eqs = eqs
     # do not need to topsort because we didn't modify the order
-    @set! sys.observed = obs
-    @set! sys.continuous_events = new_callbacks
-    @set! sys.ps = new_ps
+    @reset sys.observed = obs
+    @reset sys.continuous_events = new_callbacks
+    @reset sys.ps = new_ps
     return sys
 end
diff --git a/src/systems/imperative_affect.jl b/src/systems/imperative_affect.jl
index 2f489913ad..20db2d9dc2 100644
--- a/src/systems/imperative_affect.jl
+++ b/src/systems/imperative_affect.jl
@@ -21,10 +21,10 @@ then the NamedTuple `(;x=2)` will be passed as `observed` to the affect function
 The NamedTuple returned from `f` includes the values to be written back to the system after `f` returns. For example, if we want to update the value of `x` to be the result of `x + y` we could write
 
     ImperativeAffect(observed=(; x_plus_y = x + y), modified=(; x)) do m, o
-        @set! m.x = o.x_plus_y
+        @reset m.x = o.x_plus_y
     end
 
-Where we use Setfield to copy the tuple `m` with a new value for `x`, then return the modified value of `m`. All values updated by the tuple must have names originally declared in
+Where we use Accessors to copy the tuple `m` with a new value for `x`, then return the modified value of `m`. All values updated by the tuple must have names originally declared in
 `modified`; a runtime error will be produced if a value is written that does not appear in `modified`. The user can dynamically decide not to write a value back by not including it
 in the returned tuple, in which case the associated field will not be updated.
 """
diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl
index 4fc2f18bfd..de93fc0a0a 100644
--- a/src/systems/index_cache.jl
+++ b/src/systems/index_cache.jl
@@ -616,6 +616,6 @@ function subset_unknowns_observed(
         end
     end
     ic = @set ic.unknown_idx = unknown_idx
-    @set! ic.observed_syms_to_timeseries = observed_syms_to_timeseries
+    @reset ic.observed_syms_to_timeseries = observed_syms_to_timeseries
     return ic
 end
diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl
index 5402659b1c..2da15920d1 100644
--- a/src/systems/model_parsing.jl
+++ b/src/systems/model_parsing.jl
@@ -122,15 +122,15 @@ function _model_macro(mod, name, expr, isconnector)
     end
 
     isconnector && push!(exprs.args,
-        :($Setfield.@set!(var"#___sys___".connector_type=$connector_type(var"#___sys___"))))
+        :($Accessors.@reset(var"#___sys___".connector_type=$connector_type(var"#___sys___"))))
 
     !isempty(c_evts) && push!(exprs.args,
-        :($Setfield.@set!(var"#___sys___".continuous_events=$SymbolicContinuousCallback.([
+        :($Accessors.@reset(var"#___sys___".continuous_events=$SymbolicContinuousCallback.([
             $(c_evts...)
         ]))))
 
     !isempty(d_evts) && push!(exprs.args,
-        :($Setfield.@set!(var"#___sys___".discrete_events=$SymbolicDiscreteCallback.([
+        :($Accessors.@reset(var"#___sys___".discrete_events=$SymbolicDiscreteCallback.([
             $(d_evts...)
         ]))))
 
diff --git a/src/systems/nonlinear/homotopy_continuation.jl b/src/systems/nonlinear/homotopy_continuation.jl
index 03aeed1edf..90279202e3 100644
--- a/src/systems/nonlinear/homotopy_continuation.jl
+++ b/src/systems/nonlinear/homotopy_continuation.jl
@@ -479,10 +479,10 @@ function transform_system(sys::NonlinearSystem, transformation::PolynomialTransf
     end
 
     sys2 = @set sys.eqs = eqs2
-    @set! sys2.unknowns = new_dvs
+    @reset sys2.unknowns = new_dvs
     # remove observed equations to avoid adding them in codegen
-    @set! sys2.observed = Equation[]
-    @set! sys2.substitutions = nothing
+    @reset sys2.observed = Equation[]
+    @reset sys2.substitutions = nothing
     return PolynomialTransformationResult(sys2, denoms)
 end
 
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index b2abac5184..1ae46b3129 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -595,9 +595,9 @@ function SCCNonlinearFunction{iip}(
     subsys = NonlinearSystem(_eqs, _dvs, ps; observed = _obs,
         parameter_dependencies = parameter_dependencies(sys), name = nameof(sys))
     if get_index_cache(sys) !== nothing
-        @set! subsys.index_cache = subset_unknowns_observed(
+        @reset subsys.index_cache = subset_unknowns_observed(
             get_index_cache(sys), sys, _dvs, getproperty.(_obs, (:lhs,)))
-        @set! subsys.complete = true
+        @reset subsys.complete = true
     end
 
     return NonlinearFunction{iip}(f; sys = subsys)
@@ -713,8 +713,8 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::NonlinearSystem, u0map,
 
     new_dvs = dvs[reduce(vcat, var_sccs)]
     new_eqs = eqs[reduce(vcat, eq_sccs)]
-    @set! sys.unknowns = new_dvs
-    @set! sys.eqs = new_eqs
+    @reset sys.unknowns = new_dvs
+    @reset sys.eqs = new_eqs
     sys = complete(sys)
     return SCCNonlinearProblem(subprobs, explicitfuns, p, true; sys)
 end
diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
index 0398c892eb..a6fc9bec7e 100644
--- a/src/systems/optimization/optimizationsystem.jl
+++ b/src/systems/optimization/optimizationsystem.jl
@@ -761,11 +761,11 @@ function structural_simplify(sys::OptimizationSystem; split = true, kwargs...)
         cons_simplified[i] = fixpoint_sub(eq, subs)
     end
     newsts = setdiff(unknowns(sys), keys(subs))
-    @set! sys.constraints = cons_simplified
-    @set! sys.observed = [observed(sys); obs]
+    @reset sys.constraints = cons_simplified
+    @reset sys.observed = [observed(sys); obs]
     neweqs = fixpoint_sub.(equations(sys), (subs,))
-    @set! sys.op = length(neweqs) == 1 ? first(neweqs) : neweqs
-    @set! sys.unknowns = newsts
+    @reset sys.op = length(neweqs) == 1 ? first(neweqs) : neweqs
+    @reset sys.unknowns = newsts
     sys = complete(sys; split)
     return sys
 end
diff --git a/src/systems/parameter_buffer.jl b/src/systems/parameter_buffer.jl
index 7d64054acc..a261debcb5 100644
--- a/src/systems/parameter_buffer.jl
+++ b/src/systems/parameter_buffer.jl
@@ -294,7 +294,7 @@ function SciMLStructures.canonicalize(::SciMLStructures.Tunable, p::MTKParameter
 end
 
 function SciMLStructures.replace(::SciMLStructures.Tunable, p::MTKParameters, newvals)
-    @set! p.tunable = newvals
+    @reset p.tunable = newvals
     return p
 end
 
@@ -318,7 +318,7 @@ for (Portion, field, recurse) in [(SciMLStructures.Discrete, :discrete, 1)
     end
 
     @eval function SciMLStructures.replace(::$Portion, p::MTKParameters, newvals)
-        @set! p.$field = split_into_buffers(newvals, p.$field, Val($recurse))
+        @reset p.$field = split_into_buffers(newvals, p.$field, Val($recurse))
         p
     end
 
@@ -505,9 +505,9 @@ function SymbolicIndexingInterface.remake_buffer(indp, oldbuf::MTKParameters, id
 end
 function _remake_buffer(indp, oldbuf::MTKParameters, idxs, vals; validate = true)
     newbuf = @set oldbuf.tunable = similar(oldbuf.tunable, Any)
-    @set! newbuf.discrete = Tuple(similar(buf, Any) for buf in newbuf.discrete)
-    @set! newbuf.constant = Tuple(similar(buf, Any) for buf in newbuf.constant)
-    @set! newbuf.nonnumeric = Tuple(similar(buf, Any) for buf in newbuf.nonnumeric)
+    @reset newbuf.discrete = Tuple(similar(buf, Any) for buf in newbuf.discrete)
+    @reset newbuf.constant = Tuple(similar(buf, Any) for buf in newbuf.constant)
+    @reset newbuf.nonnumeric = Tuple(similar(buf, Any) for buf in newbuf.nonnumeric)
 
     function handle_parameter(ic, sym, idx, val)
         if sym === nothing
@@ -572,13 +572,13 @@ function _remake_buffer(indp, oldbuf::MTKParameters, idxs, vals; validate = true
         end
     end
 
-    @set! newbuf.tunable = narrow_buffer_type_and_fallback_undefs(
+    @reset newbuf.tunable = narrow_buffer_type_and_fallback_undefs(
         oldbuf.tunable, newbuf.tunable)
-    @set! newbuf.discrete = narrow_buffer_type_and_fallback_undefs.(
+    @reset newbuf.discrete = narrow_buffer_type_and_fallback_undefs.(
         oldbuf.discrete, newbuf.discrete)
-    @set! newbuf.constant = narrow_buffer_type_and_fallback_undefs.(
+    @reset newbuf.constant = narrow_buffer_type_and_fallback_undefs.(
         oldbuf.constant, newbuf.constant)
-    @set! newbuf.nonnumeric = narrow_buffer_type_and_fallback_undefs.(
+    @reset newbuf.nonnumeric = narrow_buffer_type_and_fallback_undefs.(
         oldbuf.nonnumeric, newbuf.nonnumeric)
     return newbuf
 end
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index 47acd81a82..c06ea021f6 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -50,7 +50,7 @@ function structural_simplify(
         newsys = pass(newsys)
     end
     if newsys isa ODESystem || has_parent(newsys)
-        @set! newsys.parent = complete(sys; split, flatten = false)
+        @reset newsys.parent = complete(sys; split, flatten = false)
     end
     newsys = complete(newsys; split)
     if has_defaults(newsys) && (defs = get_defaults(newsys)) !== nothing
@@ -125,8 +125,8 @@ function __structural_simplify(sys::AbstractSystem, io = nothing; simplify = fal
         end
         g = Matrix(sparse(Is, Js, vals))
         sys = state.sys
-        @set! sys.eqs = new_eqs
-        @set! sys.unknowns = [v
+        @reset sys.eqs = new_eqs
+        @reset sys.unknowns = [v
                               for (i, v) in enumerate(fullvars)
                               if !iszero(new_idxs[i]) &&
                                  invview(var_to_diff)[i] === nothing]
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index 1c0ad8b8ae..43be117f25 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -14,7 +14,7 @@ using ..BipartiteGraphs
 import ..BipartiteGraphs: invview, complete
 using Graphs
 using UnPack
-using Setfield
+using Accessors
 using SparseArrays
 
 function quick_cancel_expr(expr)
@@ -206,8 +206,8 @@ end
 TransformationState(sys::AbstractSystem) = TearingState(sys)
 function system_subset(ts::TearingState, ieqs::Vector{Int})
     eqs = equations(ts)
-    @set! ts.sys.eqs = eqs[ieqs]
-    @set! ts.structure = system_subset(ts.structure, ieqs)
+    @reset ts.sys.eqs = eqs[ieqs]
+    @reset ts.structure = system_subset(ts.structure, ieqs)
     ts
 end
 
@@ -222,8 +222,8 @@ function system_subset(structure::SystemStructure, ieqs::Vector{Int})
         push!(fadj, ivars)
         eq_to_diff[j] = structure.eq_to_diff[eq_i]
     end
-    @set! structure.graph = complete(BipartiteGraph(ne, fadj, ndsts(graph)))
-    @set! structure.eq_to_diff = eq_to_diff
+    @reset structure.graph = complete(BipartiteGraph(ne, fadj, ndsts(graph)))
+    @reset structure.eq_to_diff = eq_to_diff
     structure
 end
 
@@ -424,7 +424,7 @@ function TearingState(sys; quick_cancel = false, check = true)
         add_edge!(graph, ie, jv)
     end
 
-    @set! sys.eqs = eqs
+    @reset sys.eqs = eqs
 
     eq_to_diff = DiffGraph(nsrcs(graph))
 
@@ -475,8 +475,8 @@ function shift_discrete_system(ts::TearingState)
         eqs[i] = StructuralTransformations.simplify_shifts(fast_substitute(
             eqs[i], discmap; operator = Union{Sample, Hold}))
     end
-    @set! ts.sys.eqs = eqs
-    @set! ts.fullvars = fullvars
+    @reset ts.sys.eqs = eqs
+    @reset ts.fullvars = fullvars
     return ts
 end
 
@@ -655,16 +655,16 @@ function structural_simplify!(state::TearingState, io = nothing; simplify = fals
                 append!(appended_parameters, inputs[i], unknowns(ss))
                 discrete_subsystems[i] = ss
             end
-            @set! sys.discrete_subsystems = discrete_subsystems, inputs, continuous_id,
+            @reset sys.discrete_subsystems = discrete_subsystems, inputs, continuous_id,
             id_to_clock
-            @set! sys.ps = appended_parameters
-            @set! sys.defaults = merge(ModelingToolkit.defaults(sys),
+            @reset sys.ps = appended_parameters
+            @reset sys.defaults = merge(ModelingToolkit.defaults(sys),
                 Dict(v => 0.0 for v in Iterators.flatten(inputs)))
         end
         ps = [sym isa CallWithMetadata ? sym :
               setmetadata(sym, VariableTimeDomain, get(time_domains, sym, Continuous))
               for sym in get_ps(sys)]
-        @set! sys.ps = ps
+        @reset sys.ps = ps
     else
         sys, input_idxs = _structural_simplify!(state, io; simplify, check_consistency,
             fully_determined, kwargs...)
@@ -712,7 +712,7 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
             sys, state; simplify, mm, check_consistency, kwargs...)
     end
     fullunknowns = [map(eq -> eq.lhs, observed(sys)); unknowns(sys)]
-    @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns)
+    @reset sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns)
 
     ModelingToolkit.invalidate_cache!(sys), input_idxs
 end
diff --git a/test/clock.jl b/test/clock.jl
index 91aaa8248e..1eefa5dce5 100644
--- a/test/clock.jl
+++ b/test/clock.jl
@@ -1,4 +1,4 @@
-using ModelingToolkit, Test, Setfield, OrdinaryDiffEq, DiffEqCallbacks
+using ModelingToolkit, Test, Accessors, OrdinaryDiffEq, DiffEqCallbacks
 using ModelingToolkit: Continuous
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl
index ccf0a17a40..de1c859598 100644
--- a/test/symbolic_events.jl
+++ b/test/symbolic_events.jl
@@ -8,7 +8,7 @@ using ModelingToolkit: SymbolicContinuousCallback,
 using StableRNGs
 import SciMLBase
 using SymbolicIndexingInterface
-using Setfield
+using Accessors
 rng = StableRNG(12345)
 
 @variables x(t) = 0
@@ -1093,12 +1093,12 @@ end
     furnace_off = ModelingToolkit.SymbolicContinuousCallback(
         [temp ~ furnace_off_threshold],
         ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c
-            @set! x.furnace_on = false
+            @reset x.furnace_on = false
         end)
     furnace_enable = ModelingToolkit.SymbolicContinuousCallback(
         [temp ~ furnace_on_threshold],
         ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c
-            @set! x.furnace_on = true
+            @reset x.furnace_on = true
         end)
     @named sys = ODESystem(
         eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable])
@@ -1110,15 +1110,15 @@ end
     furnace_off = ModelingToolkit.SymbolicContinuousCallback(
         [temp ~ furnace_off_threshold],
         ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i
-            @set! x.furnace_on = false
+            @reset x.furnace_on = false
         end; initialize = ModelingToolkit.ImperativeAffect(modified = (;
             temp)) do x, o, c, i
-            @set! x.temp = 0.2
+            @reset x.temp = 0.2
         end)
     furnace_enable = ModelingToolkit.SymbolicContinuousCallback(
         [temp ~ furnace_on_threshold],
         ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i
-            @set! x.furnace_on = true
+            @reset x.furnace_on = true
         end)
     @named sys = ODESystem(
         eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable])
@@ -1140,7 +1140,7 @@ end
         [temp ~ furnace_off_threshold],
         ModelingToolkit.ImperativeAffect(
             modified = (; furnace_on), observed = (; furnace_on)) do x, o, c, i
-            @set! x.furnace_on = false
+            @reset x.furnace_on = false
         end)
     @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off])
     ss = structural_simplify(sys)
@@ -1156,7 +1156,7 @@ end
         [temp ~ furnace_off_threshold],
         ModelingToolkit.ImperativeAffect(
             modified = (; furnace_on, tempsq), observed = (; furnace_on)) do x, o, c, i
-            @set! x.furnace_on = false
+            @reset x.furnace_on = false
         end)
     @named sys = ODESystem(
         eqs, t, [temp, tempsq], params; continuous_events = [furnace_off])
@@ -1169,7 +1169,7 @@ end
         [temp ~ furnace_off_threshold],
         ModelingToolkit.ImperativeAffect(modified = (; furnace_on),
             observed = (; furnace_on, not_actually_here)) do x, o, c, i
-            @set! x.furnace_on = false
+            @reset x.furnace_on = false
         end)
     @named sys = ODESystem(
         eqs, t, [temp, tempsq], params; continuous_events = [furnace_off])
@@ -1213,34 +1213,34 @@ end
     end
     qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0],
         ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, c, i
-            @set! x.hA = x.qA
-            @set! x.hB = o.qB
-            @set! x.qA = 1
-            @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
+            @reset x.hA = x.qA
+            @reset x.hB = o.qB
+            @reset x.qA = 1
+            @reset x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
             x
         end,
         affect_neg = ModelingToolkit.ImperativeAffect(
             (; qA, hA, hB, cnt), (; qB)) do x, o, c, i
-            @set! x.hA = x.qA
-            @set! x.hB = o.qB
-            @set! x.qA = 0
-            @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
+            @reset x.hA = x.qA
+            @reset x.hB = o.qB
+            @reset x.qA = 0
+            @reset x.cnt += decoder(x.hA, x.hB, x.qA, o.qB)
             x
         end; rootfind = SciMLBase.RightRootFind)
     qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0],
         ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA)) do x, o, c, i
-            @set! x.hA = o.qA
-            @set! x.hB = x.qB
-            @set! x.qB = 1
-            @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB)
+            @reset x.hA = o.qA
+            @reset x.hB = x.qB
+            @reset x.qB = 1
+            @reset x.cnt += decoder(x.hA, x.hB, o.qA, x.qB)
             x
         end,
         affect_neg = ModelingToolkit.ImperativeAffect(
             (; qB, hA, hB, cnt), (; qA)) do x, o, c, i
-            @set! x.hA = o.qA
-            @set! x.hB = x.qB
-            @set! x.qB = 0
-            @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB)
+            @reset x.hA = o.qA
+            @reset x.hB = x.qB
+            @reset x.qB = 0
+            @reset x.cnt += decoder(x.hA, x.hB, o.qA, x.qB)
             x
         end; rootfind = SciMLBase.RightRootFind)
     @named sys = ODESystem(