diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 1987cb9a60..eaadef9be9 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -717,6 +717,12 @@ function add_initialization_parameters(sys::AbstractSystem)
             push!(all_initialvars, x)
         end
     end
+
+    # add derivatives of all variables for steady-state initial conditions
+    if is_time_dependent(sys) && !(sys isa AbstractDiscreteSystem)
+        D = Differential(get_iv(sys))
+        union!(all_initialvars, [D(v) for v in all_initialvars if iscall(v)])
+    end
     for eq in parameter_dependencies(sys)
         is_variable_floatingpoint(eq.lhs) || continue
         push!(all_initialvars, eq.lhs)
@@ -1924,7 +1930,9 @@ function toexpr(sys::AbstractSystem)
     end
 
     eqs_name = push_eqs!(stmt, full_equations(sys), var2name)
-    defs_name = push_defaults!(stmt, defaults(sys), var2name)
+    filtered_defs = filter(
+        kvp -> !(iscall(kvp[1]) && operation(kvp[1]) isa Initial), defaults(sys))
+    defs_name = push_defaults!(stmt, filtered_defs, var2name)
     obs_name = push_eqs!(stmt, obs, var2name)
 
     if sys isa ODESystem
diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl
index d0ae0d174f..1db153781a 100644
--- a/src/systems/problem_utils.jl
+++ b/src/systems/problem_utils.jl
@@ -623,12 +623,12 @@ function build_operating_point!(sys::AbstractSystem,
     end
 
     for k in keys(u0map)
-        v = fixpoint_sub(u0map[k], neithermap)
+        v = fixpoint_sub(u0map[k], neithermap; operator = Symbolics.Operator)
         isequal(k, v) && continue
         u0map[k] = v
     end
     for k in keys(pmap)
-        v = fixpoint_sub(pmap[k], neithermap)
+        v = fixpoint_sub(pmap[k], neithermap; operator = Symbolics.Operator)
         isequal(k, v) && continue
         pmap[k] = v
     end
diff --git a/test/initial_values.jl b/test/initial_values.jl
index dbe08962cd..f63333cb2f 100644
--- a/test/initial_values.jl
+++ b/test/initial_values.jl
@@ -236,3 +236,19 @@ end
     prob2 = remake(prob; p = p_new)
     @test prob2.ps[interp] == spline2
 end
+
+@testset "Issue#3523: don't substitute inside initial in `build_operating_point!`" begin
+    @variables (X(t))[1:2]
+    @parameters p[1:2]
+    eqs = [
+        0 ~ p[1] - X[1],
+        0 ~ p[2] - X[2]
+    ]
+    @named nlsys = NonlinearSystem(eqs)
+    nlsys = complete(nlsys)
+
+    # Creates the `NonlinearProblem`.
+    u0 = [X => [1.0, 2.0]]
+    ps = [p => [4.0, 5.0]]
+    @test_nowarn NonlinearProblem(nlsys, u0, ps)
+end
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 94a1f4c14f..0f40d4eaf1 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -406,30 +406,35 @@ sol = solve(prob, Tsit5())
 @test sol.u[1] == [1.0]
 
 # Steady state initialization
+@testset "Steady state initialization" begin
+    @parameters σ ρ β
+    @variables x(t) y(t) z(t)
 
-@parameters σ ρ β
-@variables x(t) y(t) z(t)
+    eqs = [D(D(x)) ~ σ * (y - x),
+        D(y) ~ x * (ρ - z) - y,
+        D(z) ~ x * y - β * z]
 
-eqs = [D(D(x)) ~ σ * (y - x),
-    D(y) ~ x * (ρ - z) - y,
-    D(z) ~ x * y - β * z]
+    @named sys = ODESystem(eqs, t)
+    sys = structural_simplify(sys)
 
-@named sys = ODESystem(eqs, t)
-sys = structural_simplify(sys)
+    u0 = [D(x) => 2.0,
+        x => 1.0,
+        D(y) => 0.0,
+        z => 0.0]
 
-u0 = [D(x) => 2.0,
-    x => 1.0,
-    D(y) => 0.0,
-    z => 0.0]
+    p = [σ => 28.0,
+        ρ => 10.0,
+        β => 8 / 3]
 
-p = [σ => 28.0,
-    ρ => 10.0,
-    β => 8 / 3]
+    tspan = (0.0, 0.2)
+    prob_mtk = ODEProblem(sys, u0, tspan, p)
+    sol = solve(prob_mtk, Tsit5())
+    @test sol[x * (ρ - z) - y][1] == 0.0
 
-tspan = (0.0, 0.2)
-prob_mtk = ODEProblem(sys, u0, tspan, p)
-sol = solve(prob_mtk, Tsit5())
-@test sol[x * (ρ - z) - y][1] == 0.0
+    prob_mtk.ps[Initial(D(y))] = 1.0
+    sol = solve(prob_mtk, Tsit5())
+    @test sol[x * (ρ - z) - y][1] == 1.0
+end
 
 @variables x(t) y(t) z(t)
 @parameters α=1.5 β=1.0 γ=3.0 δ=1.0
diff --git a/test/modelingtoolkitize.jl b/test/modelingtoolkitize.jl
index 28ce7b5e01..db99cc91a4 100644
--- a/test/modelingtoolkitize.jl
+++ b/test/modelingtoolkitize.jl
@@ -460,7 +460,7 @@ sys = modelingtoolkitize(prob)
     p = [10.0, 28.0, 2.66]
     sprob = SDEProblem(sdef!, sdeg!, u0, tspan, p)
     sys = complete(modelingtoolkitize(sprob))
-    @test length(ModelingToolkit.defaults(sys)) == 2length(u0) + length(p)
+    @test length(ModelingToolkit.defaults(sys)) == 3length(u0) + length(p)
     sprob2 = SDEProblem(sys, [], tspan)
 
     truevals = similar(u0)
diff --git a/test/symbolic_indexing_interface.jl b/test/symbolic_indexing_interface.jl
index 5979cd71ac..8b3da5fd72 100644
--- a/test/symbolic_indexing_interface.jl
+++ b/test/symbolic_indexing_interface.jl
@@ -24,7 +24,8 @@ using SciMLStructures: Tunable
         (odesys,), [x, y, t, ParameterIndex(Tunable(), 1), :x, :y]) ==
           [nothing, nothing, nothing, ParameterIndex(Tunable(), 1), nothing, nothing]
     @test isequal(
-        Set(parameter_symbols(odesys)), Set([a, b, Initial(x), Initial(y), Initial(xy)]))
+        Set(parameter_symbols(odesys)), Set([a, b, Initial(x), Initial(y), Initial(xy),
+            Initial(D(x)), Initial(D(y)), Initial(D(xy))]))
     @test all(is_independent_variable.((odesys,), [t, :t]))
     @test all(.!is_independent_variable.((odesys,), [x, y, a, :x, :y, :a]))
     @test isequal(independent_variable_symbols(odesys), [t])