diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 8b4901a19c..7379e23bc1 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -1,4 +1,5 @@
 using ModelingToolkit, OrdinaryDiffEq, NonlinearSolve, Test
+using SymbolicIndexingInterface
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @parameters g
@@ -493,27 +494,69 @@ end
         @test init(prob, Tsit5()).ps[sym] ≈ val
         @test solve(prob, Tsit5()).ps[sym] ≈ val
     end
+    function test_initializesystem(sys, u0map, pmap, p, equation)
+        isys = ModelingToolkit.generate_initializesystem(sys; u0map, pmap, guesses = ModelingToolkit.guesses(sys))
+        @test is_variable(isys, p)
+        @test equation in equations(isys) || (0 ~ -equation.rhs) in equations(isys)
+    end
     @variables x(t) y(t)
-    @parameters p
-    _p = ModelingToolkit.setdefault(p, missing)
-    _p = ModelingToolkit.setguess(_p, 0.0)
-    @mtkbuild sys = ODESystem([D(x) ~ x, _p ~ x + y], t)
-    prob = ODEProblem(sys, [x => 1.0, y => 1.0], (0.0, 1.0))
-    test_parameter(prob, _p, 2.0)
+    @parameters p q
+    u0map = Dict(x => 1.0, y => 1.0)
+    pmap = Dict()
+    pmap[q] = 1.0
+    # `missing` default, equation from ODEProblem
+    @mtkbuild sys = ODESystem([D(x) ~ x * q, D(y) ~ y * p], t; defaults = [p => missing])
+    pmap[p] = 2q
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
+    test_parameter(prob, p, 2.0)
+    # `missing` default, provided guess
     @mtkbuild sys = ODESystem(
         [D(x) ~ x, p ~ x + y], t; defaults = [p => missing], guesses = [p => 0.0])
-    prob = ODEProblem(sys, [x => 1.0, y => 1.0], (0.0, 1.0))
+    prob = ODEProblem(sys, u0map, (0.0, 1.0))
     test_parameter(prob, p, 2.0)
-    @mtkbuild sys = ODESystem([D(x) ~ x, p ~ x + y], t; guesses = [p => 0.0])
-    prob = ODEProblem(sys, [x => 1.0, y => 1.0], (0.0, 1.0), [p => missing])
+    test_initializesystem(sys, u0map, pmap, p, 0 ~ p - x - y)
+
+    # `missing` to ODEProblem, equation from default
+    @mtkbuild sys = ODESystem([D(x) ~ x * q, D(y) ~ y * p], t; defaults = [p => 2q])
+    pmap[p] = missing
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
+    test_parameter(prob, p, 2.0)
+    test_initializesystem(sys, u0map, pmap, p, 0 ~ 2q - p)
+    test_parameter(prob2, p, 2.0)
+    # `missing` to ODEProblem, provided guess
+    @mtkbuild sys = ODESystem(
+        [D(x) ~ x, p ~ x + y], t; guesses = [p => 0.0])
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
     test_parameter(prob, p, 2.0)
+    test_initializesystem(sys, u0map, pmap, p, 0 ~ x + y - p)
 
-    @mtkbuild sys = ODESystem([D(x) ~ x, p ~ x + y], t; guesses = [p => 1.0])
+    # No `missing`, default and guess
+    @mtkbuild sys = ODESystem([D(x) ~ x * q, D(y) ~ y * p], t; defaults = [p => 2q], guesses = [p => 0.0])
+    delete!(pmap, p)
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), pmap)
+    test_parameter(prob, p, 2.0)
+    test_initializesystem(sys, u0map, pmap, p, 0 ~ 2q - p)
+
+    # Should not be solved for:
+
+    # ODEProblem value with guess, no `missing`
+    @mtkbuild sys = ODESystem([D(x) ~ x * q, D(y) ~ y * p], t; guesses = [p => 0.0])
+    _pmap = merge(pmap, Dict(p => 3q))
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), _pmap)
+    @test prob.ps[p] ≈ 3.0
+    @test prob.f.initializeprob === nothing
+    # Default overridden by ODEProblem, guess provided
+    @mtkbuild sys = ODESystem([D(x) ~ q * x, D(y) ~ y * p], t; defaults = [p => 2q], guesses = [p => 1.0])
+    prob = ODEProblem(sys, u0map, (0.0, 1.0), _pmap)
+    @test prob.ps[p] ≈ 3.0
+    @test prob.f.initializeprob === nothing
+
+    @mtkbuild sys = ODESystem([D(x) ~ x, p ~ x + y], t; guesses = [p => 0.0])
     @test_throws ModelingToolkit.MissingParametersError ODEProblem(
         sys, [x => 1.0, y => 1.0], (0.0, 1.0))
     @mtkbuild sys = ODESystem([D(x) ~ x, p ~ x + y], t)
-    @test_throws ["Invalid setup", "parameter p", "guess"] ODEProblem(
-        sys, [x => 1.0, y => 1.0], (0.0, 1.0), [p => missing])
+    @test_throws ["Invalid setup", "parameter p", "guess"] ModelingToolkit.generate_initializesystem(
+        sys; u0map = Dict(x => 1.0, y => 1.0), pmap = Dict(p => missing), check_defguess = true)
 
     @testset "Null system" begin
         @variables x(t) y(t) s(t)
@@ -526,14 +569,15 @@ end
     using ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica: Fixed, Mass,
                                                                            Spring, Force,
                                                                            Damper
+    using ModelingToolkitStandardLibrary.Mechanical: TranslationalModelica as TM
     using ModelingToolkitStandardLibrary.Blocks: Constant
 
-    @named mass = Mass(; m = 1.0, s = 1.0, v = 0.0, a = 0.0)
+    @named mass = TM.Mass(; m = 1.0, s = 1.0, v = 0.0, a = 0.0)
     @named fixed = Fixed(; s0 = 0.0)
-    @named spring = Spring(; c = 2.0)
+    @named spring = Spring(; c = 2.0, s_rel0 = nothing)
     @named gravity = Force()
     @named constant = Constant(; k = 9.81)
-    @named damper = Damper(; d = 0.1)
+    @named damper = TM.Damper(; d = 0.1)
     @mtkbuild sys = ODESystem(
         [connect(fixed.flange, spring.flange_a), connect(spring.flange_b, mass.flange_a),
             connect(mass.flange_a, gravity.flange), connect(constant.output, gravity.f),