diff --git a/Project.toml b/Project.toml
index 656fef6173..67401c601a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -131,7 +131,7 @@ SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
 SymbolicIndexingInterface = "0.3.31"
 SymbolicUtils = "3.7"
-Symbolics = "6.15.2"
+Symbolics = "6.15.4"
 URIs = "1"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"
diff --git a/docs/src/examples/perturbation.md b/docs/src/examples/perturbation.md
index b2ee4bef32..cd1e843f55 100644
--- a/docs/src/examples/perturbation.md
+++ b/docs/src/examples/perturbation.md
@@ -5,20 +5,24 @@ In the [Mixed Symbolic-Numeric Perturbation Theory tutorial](https://symbolics.j
 ## Free fall in a varying gravitational field
 
 Our first ODE example is a well-known physics problem: what is the altitude $x(t)$ of an object (say, a ball or a rocket) thrown vertically with initial velocity $ẋ(0)$ from the surface of a planet with mass $M$ and radius $R$? According to Newton's second law and law of gravity, it is the solution of the ODE
+
 ```math
 ẍ = -\frac{GM}{(R+x)^2} = -\frac{GM}{R^2} \frac{1}{\left(1+ϵ\frac{x}{R}\right)^2}.
 ```
+
 In the last equality, we introduced a perturbative expansion parameter $ϵ$. When $ϵ=1$, we recover the original problem. When $ϵ=0$, the problem reduces to the trivial problem $ẍ = -g$ with constant gravitational acceleration $g = GM/R^2$ and solution $x(t) = x(0) + ẋ(0) t - \frac{1}{2} g t^2$. This is a good setup for perturbation theory.
 
 To make the problem dimensionless, we redefine $x \leftarrow x / R$ and $t \leftarrow t / \sqrt{R^3/GM}$. Then the ODE becomes
+
 ```@example perturbation
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 @variables ϵ x(t)
-eq = D(D(x)) ~ -(1 + ϵ*x)^(-2)
+eq = D(D(x)) ~ -(1 + ϵ * x)^(-2)
 ```
 
 Next, expand $x(t)$ in a series up to second order in $ϵ$:
+
 ```@example perturbation
 using Symbolics
 @variables y(t)[0:2] # coefficients
@@ -26,25 +30,32 @@ x_series = series(y, ϵ)
 ```
 
 Insert this into the equation and collect perturbed equations to each order:
+
 ```@example perturbation
 eq_pert = substitute(eq, x => x_series)
 eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2)
 ```
+
 !!! note
+    
     The 0-th order equation can be solved analytically, but ModelingToolkit does currently not feature automatic analytical solution of ODEs, so we proceed with solving it numerically.
-These are the ODEs we want to solve. Now construct an `ODESystem`, which automatically inserts dummy derivatives for the velocities:
+    These are the ODEs we want to solve. Now construct an `ODESystem`, which automatically inserts dummy derivatives for the velocities:
+
 ```@example perturbation
 @mtkbuild sys = ODESystem(eqs_pert, t)
 ```
 
 To solve the `ODESystem`, we generate an `ODEProblem` with initial conditions $x(0) = 0$, and $ẋ(0) = 1$, and solve it:
+
 ```@example perturbation
 using DifferentialEquations
 u0 = Dict([unknowns(sys) .=> 0.0; D(y[0]) => 1.0]) # nonzero initial velocity
 prob = ODEProblem(sys, u0, (0.0, 3.0))
 sol = solve(prob)
 ```
+
 This is the solution for the coefficients in the series for $x(t)$ and their derivatives. Finally, we calculate the solution to the original problem by summing the series for different $ϵ$:
+
 ```@example perturbation
 using Plots
 p = plot()
@@ -53,6 +64,7 @@ for ϵᵢ in 0.0:0.1:1.0
 end
 p
 ```
+
 This makes sense: for larger $ϵ$, gravity weakens with altitude, and the trajectory goes higher for a fixed initial velocity.
 
 An advantage of the perturbative method is that we run the ODE solver only once and calculate trajectories for several $ϵ$ for free. Had we solved the full unperturbed ODE directly, we would need to do repeat it for every $ϵ$.
@@ -62,12 +74,15 @@ An advantage of the perturbative method is that we run the ODE solver only once
 Our second example applies perturbation theory to nonlinear oscillators -- a very important class of problems. As we will see, perturbation theory has difficulty providing a good solution to this problem, but the process is nevertheless instructive. This example closely follows chapter 7.6 of *Nonlinear Dynamics and Chaos* by Steven Strogatz.
 
 The goal is to solve the ODE
+
 ```@example perturbation
-eq = D(D(x)) + 2*ϵ*D(x) + x ~ 0
+eq = D(D(x)) + 2 * ϵ * D(x) + x ~ 0
 ```
+
 with initial conditions $x(0) = 0$ and $ẋ(0) = 1$. With $ϵ = 0$, the problem reduces to the simple linear harmonic oscillator with the exact solution $x(t) = \sin(t)$.
 
 We follow the same steps as in the previous example to construct the `ODESystem`:
+
 ```@example perturbation
 eq_pert = substitute(eq, x => x_series)
 eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2)
@@ -75,6 +90,7 @@ eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2)
 ```
 
 We solve and plot it as in the previous example, and compare the solution with $ϵ=0.1$ to the exact solution $x(t, ϵ) = e^{-ϵ t} \sin(\sqrt{(1-ϵ^2)}\,t) / \sqrt{1-ϵ^2}$ of the unperturbed equation:
+
 ```@example perturbation
 u0 = Dict([unknowns(sys) .=> 0.0; D(y[0]) => 1.0]) # nonzero initial velocity
 prob = ODEProblem(sys, u0, (0.0, 50.0))
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 10401d91c5..9dda345f14 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -3001,8 +3001,8 @@ By default, the resulting system inherits `sys`'s name and description.
 See also [`compose`](@ref).
 """
 function extend(sys::AbstractSystem, basesys::AbstractSystem;
-    name::Symbol = nameof(sys), description = description(sys),
-    gui_metadata = get_gui_metadata(sys))
+        name::Symbol = nameof(sys), description = description(sys),
+        gui_metadata = get_gui_metadata(sys))
     T = SciMLBase.parameterless_type(basesys)
     ivs = independent_variables(basesys)
     if !(sys isa T)
diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl
index 977f637d47..142d2f6d71 100644
--- a/src/systems/problem_utils.jl
+++ b/src/systems/problem_utils.jl
@@ -250,6 +250,25 @@ function add_parameter_dependencies!(sys::AbstractSystem, varmap::AbstractDict)
     add_observed_equations!(varmap, parameter_dependencies(sys))
 end
 
+struct UnexpectedSymbolicValueInVarmap <: Exception
+    sym::Any
+    val::Any
+end
+
+function Base.showerror(io::IO, err::UnexpectedSymbolicValueInVarmap)
+    println(io,
+        """
+        Found symbolic value $(err.val) for variable $(err.sym). You may be missing an \
+        initial condition or have cyclic initial conditions. If this is intended, pass \
+        `symbolic_u0 = true`. In case the initial conditions are not cyclic but \
+        require more substitutions to resolve, increase `substitution_limit`. To report \
+        cycles in initial conditions of unknowns/parameters, pass \
+        `warn_cyclic_dependency = true`. If the cycles are still not reported, you \
+        may need to pass a larger value for `circular_dependency_max_cycle_length` \
+        or `circular_dependency_max_cycles`.
+        """)
+end
+
 """
     $(TYPEDSIGNATURES)
 
@@ -278,6 +297,12 @@ function better_varmap_to_vars(varmap::AbstractDict, vars::Vector;
         isempty(missing_vars) || throw(MissingVariablesError(missing_vars))
     end
     vals = map(x -> varmap[x], vars)
+    if !allow_symbolic
+        for (sym, val) in zip(vars, vals)
+            symbolic_type(val) == NotSymbolic() && continue
+            throw(UnexpectedSymbolicValueInVarmap(sym, val))
+        end
+    end
 
     if container_type <: Union{AbstractDict, Tuple, Nothing, SciMLBase.NullParameters}
         container_type = Array
@@ -298,6 +323,57 @@ function better_varmap_to_vars(varmap::AbstractDict, vars::Vector;
     end
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Check if any of the substitution rules in `varmap` lead to cycles involving
+variables in `vars`. Return a vector of vectors containing all the variables
+in each cycle.
+
+Keyword arguments:
+- `max_cycle_length`: The maximum length (number of variables) of detected cycles.
+- `max_cycles`: The maximum number of cycles to report.
+"""
+function check_substitution_cycles(
+        varmap::AbstractDict, vars; max_cycle_length = length(varmap), max_cycles = 10)
+    # ordered set so that `vars` are the first `k` in the list
+    allvars = OrderedSet{Any}(vars)
+    union!(allvars, keys(varmap))
+    allvars = collect(allvars)
+    var_to_idx = Dict(allvars .=> eachindex(allvars))
+    graph = SimpleDiGraph(length(allvars))
+
+    buffer = Set()
+    for (k, v) in varmap
+        kidx = var_to_idx[k]
+        if symbolic_type(v) != NotSymbolic()
+            vars!(buffer, v)
+            for var in buffer
+                haskey(var_to_idx, var) || continue
+                add_edge!(graph, kidx, var_to_idx[var])
+            end
+        elseif v isa AbstractArray
+            for val in v
+                vars!(buffer, val)
+            end
+            for var in buffer
+                haskey(var_to_idx, var) || continue
+                add_edge!(graph, kidx, var_to_idx[var])
+            end
+        end
+        empty!(buffer)
+    end
+
+    # detect at most 100 cycles involving at most `length(varmap)` vertices
+    cycles = Graphs.simplecycles_limited_length(graph, max_cycle_length, max_cycles)
+    # only count those which contain variables in `vars`
+    filter!(Base.Fix1(any, <=(length(vars))), cycles)
+
+    map(cycles) do cycle
+        map(Base.Fix1(getindex, allvars), cycle)
+    end
+end
+
 """
     $(TYPEDSIGNATURES)
 
@@ -305,10 +381,10 @@ Performs symbolic substitution on the values in `varmap` for the keys in `vars`,
 `varmap` itself as the set of substitution rules. If an entry in `vars` is not a key
 in `varmap`, it is ignored.
 """
-function evaluate_varmap!(varmap::AbstractDict, vars)
+function evaluate_varmap!(varmap::AbstractDict, vars; limit = 100)
     for k in vars
         haskey(varmap, k) || continue
-        varmap[k] = fixpoint_sub(varmap[k], varmap)
+        varmap[k] = fixpoint_sub(varmap[k], varmap; maxiters = limit)
     end
 end
 
@@ -407,6 +483,14 @@ Keyword arguments:
   length of `u0` vector for consistency. If `false`, do not check with equations. This is
   forwarded to `check_eqs_u0`
 - `symbolic_u0` allows the returned `u0` to be an array of symbolics.
+- `warn_cyclic_dependency`: Whether to emit a warning listing out cycles in initial
+  conditions provided for unknowns and parameters.
+- `circular_dependency_max_cycle_length`: Maximum length of cycle to check for.
+  Only applicable if `warn_cyclic_dependency == true`.
+- `circular_dependency_max_cycles`: Maximum number of cycles to check for.
+  Only applicable if `warn_cyclic_dependency == true`.
+- `substitution_limit`: The number times to substitute initial conditions into each
+  other to attempt to arrive at a numeric value.
 
 All other keyword arguments are passed as-is to `constructor`.
 """
@@ -416,7 +500,11 @@ function process_SciMLProblem(
         warn_initialize_determined = true, initialization_eqs = [],
         eval_expression = false, eval_module = @__MODULE__, fully_determined = false,
         check_initialization_units = false, tofloat = true, use_union = false,
-        u0_constructor = identity, du0map = nothing, check_length = true, symbolic_u0 = false, kwargs...)
+        u0_constructor = identity, du0map = nothing, check_length = true,
+        symbolic_u0 = false, warn_cyclic_dependency = false,
+        circular_dependency_max_cycle_length = length(all_symbols(sys)),
+        circular_dependency_max_cycles = 10,
+        substitution_limit = 100, kwargs...)
     dvs = unknowns(sys)
     ps = parameters(sys)
     iv = has_iv(sys) ? get_iv(sys) : nothing
@@ -466,7 +554,8 @@ function process_SciMLProblem(
             initializeprob = ModelingToolkit.InitializationProblem(
                 sys, t, u0map, pmap; guesses, warn_initialize_determined,
                 initialization_eqs, eval_expression, eval_module, fully_determined,
-                check_units = check_initialization_units)
+                warn_cyclic_dependency, check_units = check_initialization_units,
+                circular_dependency_max_cycle_length, circular_dependency_max_cycles)
             initializeprobmap = getu(initializeprob, unknowns(sys))
 
             punknowns = [p
@@ -503,7 +592,20 @@ function process_SciMLProblem(
     add_observed!(sys, op)
     add_parameter_dependencies!(sys, op)
 
-    evaluate_varmap!(op, dvs)
+    if warn_cyclic_dependency
+        cycles = check_substitution_cycles(
+            op, dvs; max_cycle_length = circular_dependency_max_cycle_length,
+            max_cycles = circular_dependency_max_cycles)
+        if !isempty(cycles)
+            buffer = IOBuffer()
+            for cycle in cycles
+                println(buffer, cycle)
+            end
+            msg = String(take!(buffer))
+            @warn "Cycles in unknowns:\n$msg"
+        end
+    end
+    evaluate_varmap!(op, dvs; limit = substitution_limit)
 
     u0 = better_varmap_to_vars(
         op, dvs; tofloat = true, use_union = false,
@@ -515,7 +617,20 @@ function process_SciMLProblem(
 
     check_eqs_u0(eqs, dvs, u0; check_length, kwargs...)
 
-    evaluate_varmap!(op, ps)
+    if warn_cyclic_dependency
+        cycles = check_substitution_cycles(
+            op, ps; max_cycle_length = circular_dependency_max_cycle_length,
+            max_cycles = circular_dependency_max_cycles)
+        if !isempty(cycles)
+            buffer = IOBuffer()
+            for cycle in cycles
+                println(buffer, cycle)
+            end
+            msg = String(take!(buffer))
+            @warn "Cycles in parameters:\n$msg"
+        end
+    end
+    evaluate_varmap!(op, ps; limit = substitution_limit)
     if is_split(sys)
         p = MTKParameters(sys, op)
     else
diff --git a/test/initial_values.jl b/test/initial_values.jl
index fc386242a6..a2931e22b6 100644
--- a/test/initial_values.jl
+++ b/test/initial_values.jl
@@ -149,3 +149,28 @@ end
     pmap = [c1 => 5.0, c2 => 1.0, c3 => 1.2]
     oprob = ODEProblem(osys, u0map, (0.0, 10.0), pmap)
 end
+
+@testset "Cyclic dependency checking and substitution limits" begin
+    @variables x(t) y(t)
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x, D(y) ~ y], t; initialization_eqs = [x ~ 2y + 3, y ~ 2x],
+        guesses = [x => 2y, y => 2x])
+    @test_warn ["Cycle", "unknowns", "x", "y"] try
+        ODEProblem(sys, [], (0.0, 1.0), warn_cyclic_dependency = true)
+    catch
+    end
+    @test_throws ModelingToolkit.UnexpectedSymbolicValueInVarmap ODEProblem(
+        sys, [x => 2y + 1, y => 2x], (0.0, 1.0); build_initializeprob = false)
+
+    @parameters p q
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x * p, D(y) ~ y * q], t; guesses = [p => 1.0, q => 2.0])
+    # "unknowns" because they are initialization unknowns
+    @test_warn ["Cycle", "unknowns", "p", "q"] try
+        ODEProblem(sys, [x => 1, y => 2], (0.0, 1.0),
+            [p => 2q, q => 3p]; warn_cyclic_dependency = true)
+    catch
+    end
+    @test_throws ModelingToolkit.UnexpectedSymbolicValueInVarmap ODEProblem(
+        sys, [x => 1, y => 2], (0.0, 1.0), [p => 2q, q => 3p])
+end