From 2de171c793237e84d6371f602e7ca36c804a84dc Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 21 Apr 2025 10:21:45 +0530
Subject: [PATCH 01/10] fix: copy initials to `u0` if `u0` not provided to
 `remake`

---
 src/systems/nonlinear/initializesystem.jl | 12 ++++++++++++
 test/initializationsystem.jl              | 21 +++++++++++++++++++++
 2 files changed, 33 insertions(+)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 50c0862863..8e5b59c52d 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -593,6 +593,18 @@ function SciMLBase.late_binding_update_u0_p(
         prob, sys::AbstractSystem, u0, p, t0, newu0, newp)
     supports_initialization(sys) || return newu0, newp
     u0 === missing && return newu0, (p === missing ? copy(newp) : newp)
+    # If the user passes `p` to `remake` but not `u0` and `u0` isn't empty,
+    # and if the system supports initialization (so it has initial parameters),
+    # and if the initialization solves for `u0`,
+    # THEN copy the values of `Initial`s to `newu0`.
+    if u0 === missing && newu0 !== nothing && p !== missing && supports_initialization(sys) && prob.f.initialization_data !== nothing && prob.f.initialization_data.initializeprobmap !== nothing
+        if ArrayInterface.ismutable(newu0)
+            copyto!(newu0, getu(sys, Initial.(unknowns(sys)))(newp))
+        else
+            T = StaticArrays.similar_type(newu0)
+            newu0 = T(getu(sys, Initial.(unknowns(sys)))(newp))
+        end
+    end
     # non-symbolic u0 updates initials...
     if !(eltype(u0) <: Pair)
         # if `p` is not provided or is symbolic
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 2804c70833..de7fd287ce 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -1512,3 +1512,24 @@ end
     @inferred remake(prob; u0 = 2 .* prob.u0, p = prob.p)
     @inferred solve(prob)
 end
+
+@testset "Issue#3570: `Initial`s are copied to `u0` if `u0` not provided to `remake`" begin
+    @parameters g
+    @variables x(t) [state_priority = 10] y(t) λ(t)
+    eqs = [D(D(x)) ~ λ * x
+           D(D(y)) ~ λ * y - g
+           x^2 + y^2 ~ 1]
+    @mtkbuild pend = ODESystem(eqs, t)
+
+    prob = ODEProblem(
+        pend, [x => (√2 / 2)], (0.0, 1.5), [g => 1], guesses = [λ => 1, y => √2 / 2])
+    sol = solve(prob)
+
+    setter = setsym_oop(prob, [Initial(x)])
+    (u0, p) = setter(prob, [0.8])
+
+    new_prob = remake(prob; p, initializealg = BrownFullBasicInit())
+    @test new_prob[x] ≈ 0.8
+    new_sol = solve(new_prob)
+    @test new_sol[x, 1] ≈ 0.8
+end

From eb37a14f4f407bc56b50fd8bfc3db97441014a4d Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 21 Apr 2025 15:25:28 +0530
Subject: [PATCH 02/10] fix: remove early exit in `late_binding_update_u0_p`
 when `u0 === missing`

---
 src/systems/nonlinear/initializesystem.jl | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 8e5b59c52d..a6fef3b84e 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -592,19 +592,22 @@ end
 function SciMLBase.late_binding_update_u0_p(
         prob, sys::AbstractSystem, u0, p, t0, newu0, newp)
     supports_initialization(sys) || return newu0, newp
-    u0 === missing && return newu0, (p === missing ? copy(newp) : newp)
     # If the user passes `p` to `remake` but not `u0` and `u0` isn't empty,
     # and if the system supports initialization (so it has initial parameters),
     # and if the initialization solves for `u0`,
     # THEN copy the values of `Initial`s to `newu0`.
-    if u0 === missing && newu0 !== nothing && p !== missing && supports_initialization(sys) && prob.f.initialization_data !== nothing && prob.f.initialization_data.initializeprobmap !== nothing
-        if ArrayInterface.ismutable(newu0)
-            copyto!(newu0, getu(sys, Initial.(unknowns(sys)))(newp))
-        else
-            T = StaticArrays.similar_type(newu0)
-            newu0 = T(getu(sys, Initial.(unknowns(sys)))(newp))
+    if u0 === missing
+        if newu0 !== nothing && p !== missing && supports_initialization(sys) && prob.f.initialization_data !== nothing && prob.f.initialization_data.initializeprobmap !== nothing
+            if ArrayInterface.ismutable(newu0)
+                copyto!(newu0, getu(sys, Initial.(unknowns(sys)))(newp))
+            else
+                T = StaticArrays.similar_type(newu0)
+                newu0 = T(getu(sys, Initial.(unknowns(sys)))(newp))
+            end
         end
+        return newu0, newp
     end
+
     # non-symbolic u0 updates initials...
     if !(eltype(u0) <: Pair)
         # if `p` is not provided or is symbolic

From 6aa661def8a4eb8d96afe723a5c7c0892cae0ee0 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 23 Apr 2025 20:14:53 +0530
Subject: [PATCH 03/10] refactor: store `getu` function in
 `InitializationMetadata`

---
 src/systems/nonlinear/initializesystem.jl | 16 ++++++++++++----
 src/systems/problem_utils.jl              |  9 +++++++--
 2 files changed, 19 insertions(+), 6 deletions(-)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index a6fef3b84e..332e4b1bbb 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -592,17 +592,26 @@ end
 function SciMLBase.late_binding_update_u0_p(
         prob, sys::AbstractSystem, u0, p, t0, newu0, newp)
     supports_initialization(sys) || return newu0, newp
+
+    initdata = prob.f.initialization_data
+    meta = initdata === nothing ? nothing : initdata.metadata
     # If the user passes `p` to `remake` but not `u0` and `u0` isn't empty,
     # and if the system supports initialization (so it has initial parameters),
     # and if the initialization solves for `u0`,
     # THEN copy the values of `Initial`s to `newu0`.
     if u0 === missing
-        if newu0 !== nothing && p !== missing && supports_initialization(sys) && prob.f.initialization_data !== nothing && prob.f.initialization_data.initializeprobmap !== nothing
+        if newu0 !== nothing && p !== missing && supports_initialization(sys) &&
+           initdata !== nothing && initdata.initializeprobmap !== nothing
+            getter = if meta isa InitializationMetadata
+                meta.get_initial_unknowns
+            else
+                getu(sys, Initial.(unknowns(sys)))
+            end
             if ArrayInterface.ismutable(newu0)
-                copyto!(newu0, getu(sys, Initial.(unknowns(sys)))(newp))
+                copyto!(newu0, getter(newp))
             else
                 T = StaticArrays.similar_type(newu0)
-                newu0 = T(getu(sys, Initial.(unknowns(sys)))(newp))
+                newu0 = T(getter(newp))
             end
         end
         return newu0, newp
@@ -613,7 +622,6 @@ function SciMLBase.late_binding_update_u0_p(
         # if `p` is not provided or is symbolic
         p === missing || eltype(p) <: Pair || return newu0, newp
         (newu0 === nothing || isempty(newu0)) && return newu0, newp
-        initdata = prob.f.initialization_data
         initdata === nothing && return newu0, newp
         meta = initdata.metadata
         meta isa InitializationMetadata || return newu0, newp
diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl
index 41c64b78c5..6d1a09b86a 100644
--- a/src/systems/problem_utils.jl
+++ b/src/systems/problem_utils.jl
@@ -769,7 +769,7 @@ properly.
 
 $(TYPEDFIELDS)
 """
-struct InitializationMetadata{R <: ReconstructInitializeprob, SIU}
+struct InitializationMetadata{R <: ReconstructInitializeprob, GIU, SIU}
     """
     The `u0map` used to construct the initialization.
     """
@@ -796,6 +796,11 @@ struct InitializationMetadata{R <: ReconstructInitializeprob, SIU}
     """
     oop_reconstruct_u0_p::R
     """
+    A function which takes the parameter object of the problem and returns
+    `Initial.(unknowns(sys))`.
+    """
+    get_initial_unknowns::GIU
+    """
     A function which takes the `u0` of the problem and sets
     `Initial.(unknowns(sys))`.
     """
@@ -843,7 +848,7 @@ function maybe_build_initialization_problem(
     meta = InitializationMetadata(
         u0map, pmap, guesses, Vector{Equation}(initialization_eqs),
         use_scc, ReconstructInitializeprob(sys, initializeprob.f.sys),
-        setp(sys, Initial.(unknowns(sys))))
+        getp(sys, Initial.(unknowns(sys))), setp(sys, Initial.(unknowns(sys))))
 
     if is_time_dependent(sys)
         all_init_syms = Set(all_symbols(initializeprob))

From 8b73ea21c51b900a27f832e72521f9f5bee95400 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 24 Apr 2025 11:07:47 +0530
Subject: [PATCH 04/10] refactor: copy initials to `u0` in `solve`/`init`
 instead of `remake`

---
 src/systems/nonlinear/initializesystem.jl | 55 ++++++++++++++---------
 test/initializationsystem.jl              | 25 ++++++++---
 2 files changed, 52 insertions(+), 28 deletions(-)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 332e4b1bbb..c081a02c70 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -595,27 +595,6 @@ function SciMLBase.late_binding_update_u0_p(
 
     initdata = prob.f.initialization_data
     meta = initdata === nothing ? nothing : initdata.metadata
-    # If the user passes `p` to `remake` but not `u0` and `u0` isn't empty,
-    # and if the system supports initialization (so it has initial parameters),
-    # and if the initialization solves for `u0`,
-    # THEN copy the values of `Initial`s to `newu0`.
-    if u0 === missing
-        if newu0 !== nothing && p !== missing && supports_initialization(sys) &&
-           initdata !== nothing && initdata.initializeprobmap !== nothing
-            getter = if meta isa InitializationMetadata
-                meta.get_initial_unknowns
-            else
-                getu(sys, Initial.(unknowns(sys)))
-            end
-            if ArrayInterface.ismutable(newu0)
-                copyto!(newu0, getter(newp))
-            else
-                T = StaticArrays.similar_type(newu0)
-                newu0 = T(getter(newp))
-            end
-        end
-        return newu0, newp
-    end
 
     # non-symbolic u0 updates initials...
     if !(eltype(u0) <: Pair)
@@ -669,6 +648,40 @@ function SciMLBase.late_binding_update_u0_p(
     return newu0, newp
 end
 
+function DiffEqBase.get_updated_symbolic_problem(sys::AbstractSystem, prob; kw...)
+    supports_initialization(sys) || return prob
+    initdata = prob.f.initialization_data
+    initdata === nothing && return prob
+
+    u0 = state_values(prob)
+    u0 === nothing && return prob
+
+    p = parameter_values(prob)
+    t0 = is_time_dependent(prob) ? current_time(prob) : nothing
+    meta = initdata.metadata
+
+    getter = if meta isa InitializationMetadata
+        meta.get_initial_unknowns
+    else
+        getu(sys, Initial.(unknowns(sys)))
+    end
+
+    if p isa MTKParameters
+        buffer = p.initials
+    else
+        buffer = p
+    end
+
+    u0 = DiffEqBase.promote_u0(u0, buffer, t0)
+
+    if ArrayInterface.ismutable(u0)
+        T = typeof(u0)
+    else
+        T = StaticArrays.similar_type(u0)
+    end
+    return remake(prob; u0 = T(getter(p)))
+end
+
 """
     $(TYPEDSIGNATURES)
 
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index de7fd287ce..841d29878b 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -1513,7 +1513,7 @@ end
     @inferred solve(prob)
 end
 
-@testset "Issue#3570: `Initial`s are copied to `u0` if `u0` not provided to `remake`" begin
+@testset "Issue#3570, #3552: `Initial`s are copied to `u0` during `solve`/`init`" begin
     @parameters g
     @variables x(t) [state_priority = 10] y(t) λ(t)
     eqs = [D(D(x)) ~ λ * x
@@ -1525,11 +1525,22 @@ end
         pend, [x => (√2 / 2)], (0.0, 1.5), [g => 1], guesses = [λ => 1, y => √2 / 2])
     sol = solve(prob)
 
-    setter = setsym_oop(prob, [Initial(x)])
-    (u0, p) = setter(prob, [0.8])
+    @testset "`setsym_oop`" begin
+        setter = setsym_oop(prob, [Initial(x)])
+        (u0, p) = setter(prob, [0.8])
+        new_prob = remake(prob; u0, p, initializealg = BrownFullBasicInit())
+        new_sol = solve(new_prob)
+        @test new_sol[x, 1] ≈ 0.8
+        integ = init(new_prob)
+        @test integ[x] ≈ 0.8
+    end
 
-    new_prob = remake(prob; p, initializealg = BrownFullBasicInit())
-    @test new_prob[x] ≈ 0.8
-    new_sol = solve(new_prob)
-    @test new_sol[x, 1] ≈ 0.8
+    @testset "`setsym`" begin
+        @test prob.ps[Initial(x)] ≈ √2 / 2
+        prob.ps[Initial(x)] = 0.8
+        sol = solve(prob; initializealg = BrownFullBasicInit())
+        @test sol[x, 1] ≈ 0.8
+        integ = init(prob; initializealg = BrownFullBasicInit())
+        @test integ[x] ≈ 0.8
+    end
 end

From 2bfe21fa7556a05fdb80798f019dcbe547656046 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 24 Apr 2025 15:45:22 +0530
Subject: [PATCH 05/10] build: bump DiffEqBase compat

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 35e03e2924..8a7cd147b4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -95,7 +95,7 @@ DataInterpolations = "6.4"
 DataStructures = "0.17, 0.18"
 DeepDiffs = "1"
 DelayDiffEq = "5.50"
-DiffEqBase = "6.165.1"
+DiffEqBase = "6.170.1"
 DiffEqCallbacks = "2.16, 3, 4"
 DiffEqNoiseProcess = "5"
 DiffRules = "0.1, 1.0"

From cc2f068a88dfbe3d403e5b2e626c1c9081ab69d0 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 24 Apr 2025 18:47:53 +0530
Subject: [PATCH 06/10] fix: fix type promotion in `late_binding_update_u0_p`

---
 src/systems/nonlinear/initializesystem.jl | 36 ++++++++++++-----------
 1 file changed, 19 insertions(+), 17 deletions(-)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index c081a02c70..67f99f912b 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -589,6 +589,22 @@ function SciMLBase.remake_initialization_data(
     return SciMLBase.remake_initialization_data(sys, kws, newu0, t0, newp, newu0, newp)
 end
 
+function promote_u0_p(u0, p::MTKParameters, t0)
+    u0 = DiffEqBase.promote_u0(u0, p.tunable, t0)
+    u0 = DiffEqBase.promote_u0(u0, p.initials, t0)
+
+    tunables = DiffEqBase.promote_u0(p.tunable, u0, t0)
+    initials = DiffEqBase.promote_u0(p.initials, u0, t0)
+    p = SciMLStructures.replace(SciMLStructures.Tunable(), p, tunables)
+    p = SciMLStructures.replace(SciMLStructures.Initials(), p, initials)
+
+    return u0, p
+end
+
+function promote_u0_p(u0, p::AbstractArray, t0)
+    return DiffEqBase.promote_u0(u0, p, t0), DiffEqBase.promote_u0(p, u0, t0)
+end
+
 function SciMLBase.late_binding_update_u0_p(
         prob, sys::AbstractSystem, u0, p, t0, newu0, newp)
     supports_initialization(sys) || return newu0, newp
@@ -596,6 +612,8 @@ function SciMLBase.late_binding_update_u0_p(
     initdata = prob.f.initialization_data
     meta = initdata === nothing ? nothing : initdata.metadata
 
+    newu0, newp = promote_u0_p(newu0, newp, t0)
+
     # non-symbolic u0 updates initials...
     if !(eltype(u0) <: Pair)
         # if `p` is not provided or is symbolic
@@ -605,12 +623,7 @@ function SciMLBase.late_binding_update_u0_p(
         meta = initdata.metadata
         meta isa InitializationMetadata || return newu0, newp
         newp = p === missing ? copy(newp) : newp
-        initials, repack, alias = SciMLStructures.canonicalize(
-            SciMLStructures.Initials(), newp)
-        if eltype(initials) != eltype(newu0)
-            initials = DiffEqBase.promote_u0(initials, newu0, t0)
-            newp = repack(initials)
-        end
+
         if length(newu0) != length(prob.u0)
             throw(ArgumentError("Expected `newu0` to be of same length as unknowns ($(length(prob.u0))). Got $(typeof(newu0)) of length $(length(newu0))"))
         end
@@ -619,17 +632,6 @@ function SciMLBase.late_binding_update_u0_p(
     end
 
     newp = p === missing ? copy(newp) : newp
-    newu0 = DiffEqBase.promote_u0(newu0, newp, t0)
-    tunables, repack, alias = SciMLStructures.canonicalize(SciMLStructures.Tunable(), newp)
-    if eltype(tunables) != eltype(newu0)
-        tunables = DiffEqBase.promote_u0(tunables, newu0, t0)
-        newp = repack(tunables)
-    end
-    initials, repack, alias = SciMLStructures.canonicalize(SciMLStructures.Initials(), newp)
-    if eltype(initials) != eltype(newu0)
-        initials = DiffEqBase.promote_u0(initials, newu0, t0)
-        newp = repack(initials)
-    end
 
     allsyms = all_symbols(sys)
     for (k, v) in u0

From df5afa92130bcad34e1d0ab96c98b27fc8e34fe1 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 25 Apr 2025 19:19:50 +0530
Subject: [PATCH 07/10] fix: propagate guesses for algebraic variables instead
 of `Initial(x)`

---
 src/systems/nonlinear/initializesystem.jl | 15 +++---
 src/systems/problem_utils.jl              | 56 +++++++++++++++++++++--
 2 files changed, 57 insertions(+), 14 deletions(-)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 67f99f912b..5995272be9 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -653,20 +653,16 @@ end
 function DiffEqBase.get_updated_symbolic_problem(sys::AbstractSystem, prob; kw...)
     supports_initialization(sys) || return prob
     initdata = prob.f.initialization_data
-    initdata === nothing && return prob
+    initdata isa SciMLBase.OverrideInitData || return prob
+    meta = initdata.metadata
+    meta isa InitializationMetadata || return prob
+    meta.get_updated_u0 === nothing && return prob
 
     u0 = state_values(prob)
     u0 === nothing && return prob
 
     p = parameter_values(prob)
     t0 = is_time_dependent(prob) ? current_time(prob) : nothing
-    meta = initdata.metadata
-
-    getter = if meta isa InitializationMetadata
-        meta.get_initial_unknowns
-    else
-        getu(sys, Initial.(unknowns(sys)))
-    end
 
     if p isa MTKParameters
         buffer = p.initials
@@ -681,7 +677,8 @@ function DiffEqBase.get_updated_symbolic_problem(sys::AbstractSystem, prob; kw..
     else
         T = StaticArrays.similar_type(u0)
     end
-    return remake(prob; u0 = T(getter(p)))
+
+    return remake(prob; u0 = T(meta.get_updated_u0(prob, initdata.initializeprob)))
 end
 
 """
diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl
index 6d1a09b86a..b46ef50f8e 100644
--- a/src/systems/problem_utils.jl
+++ b/src/systems/problem_utils.jl
@@ -769,7 +769,7 @@ properly.
 
 $(TYPEDFIELDS)
 """
-struct InitializationMetadata{R <: ReconstructInitializeprob, GIU, SIU}
+struct InitializationMetadata{R <: ReconstructInitializeprob, GUU, SIU}
     """
     The `u0map` used to construct the initialization.
     """
@@ -796,10 +796,9 @@ struct InitializationMetadata{R <: ReconstructInitializeprob, GIU, SIU}
     """
     oop_reconstruct_u0_p::R
     """
-    A function which takes the parameter object of the problem and returns
-    `Initial.(unknowns(sys))`.
+    A function which takes `(prob, initializeprob)` and return the `u0` to use for the problem.
     """
-    get_initial_unknowns::GIU
+    get_updated_u0::GUU
     """
     A function which takes the `u0` of the problem and sets
     `Initial.(unknowns(sys))`.
@@ -807,6 +806,48 @@ struct InitializationMetadata{R <: ReconstructInitializeprob, GIU, SIU}
     set_initial_unknowns!::SIU
 end
 
+"""
+    $(TYPEDEF)
+
+A callable struct to use as the `get_updated_u0` field of `InitializationMetadata`.
+Returns the value of `Initial.(unknowns(sys))`, except with algebraic variables replaced
+by their guess values in the initialization problem.
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
+struct GetUpdatedU0{GA, GIU}
+    """
+    Mask with length `length(unknowns(sys))` denoting indices of algebraic variables.
+    """
+    algevars::BitVector
+    """
+    Function which returns the values of algebraic variables in `initializeprob`, in the
+    order the algebraic variables occur in `unknowns(sys)`.
+    """
+    get_algevars::GA
+    """
+    Function which returns `Initial.(unknowns(sys))` as a `Vector`.
+    """
+    get_initial_unknowns::GIU
+end
+
+function GetUpdatedU0(sys::AbstractSystem, initsys::AbstractSystem)
+    algevaridxs = BitVector(is_alg_equation.(equations(sys)))
+    algevars = unknowns(sys)[algevaridxs]
+    get_algevars = getu(initsys, algevars)
+    get_initial_unknowns = getu(sys, Initial.(unknowns(sys)))
+    return GetUpdatedU0(algevaridxs, get_algevars, get_initial_unknowns)
+end
+
+function (guu::GetUpdatedU0)(prob, initprob)
+    buffer = guu.get_initial_unknowns(prob)
+    algebuf = view(buffer, guu.algevars)
+    copyto!(algebuf, guu.get_algevars(initprob))
+    return buffer
+end
+
 """
     $(TYPEDSIGNATURES)
 
@@ -845,10 +886,15 @@ function maybe_build_initialization_problem(
     end
     initializeprob = remake(initializeprob; p = initp)
 
+    get_initial_unknowns = if is_time_dependent(sys)
+        GetUpdatedU0(sys, initializeprob.f.sys)
+    else
+        nothing
+    end
     meta = InitializationMetadata(
         u0map, pmap, guesses, Vector{Equation}(initialization_eqs),
         use_scc, ReconstructInitializeprob(sys, initializeprob.f.sys),
-        getp(sys, Initial.(unknowns(sys))), setp(sys, Initial.(unknowns(sys))))
+        get_initial_unknowns, setp(sys, Initial.(unknowns(sys))))
 
     if is_time_dependent(sys)
         all_init_syms = Set(all_symbols(initializeprob))

From bc01882c588594bf51fd49b21082e0297f48f36d Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 25 Apr 2025 19:48:43 +0530
Subject: [PATCH 08/10] test: test `get_updated_symbolic_problem` copying
 guesses of algebraic variables

---
 test/initializationsystem.jl | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 841d29878b..2534fb163f 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -1513,7 +1513,7 @@ end
     @inferred solve(prob)
 end
 
-@testset "Issue#3570, #3552: `Initial`s are copied to `u0` during `solve`/`init`" begin
+@testset "Issue#3570, #3552: `Initial`s/guesses are copied to `u0` during `solve`/`init`" begin
     @parameters g
     @variables x(t) [state_priority = 10] y(t) λ(t)
     eqs = [D(D(x)) ~ λ * x
@@ -1522,9 +1522,17 @@ end
     @mtkbuild pend = ODESystem(eqs, t)
 
     prob = ODEProblem(
-        pend, [x => (√2 / 2)], (0.0, 1.5), [g => 1], guesses = [λ => 1, y => √2 / 2])
+        pend, [x => (√2 / 2), D(x) => 0.0], (0.0, 1.5),
+        [g => 1], guesses = [λ => 1, y => √2 / 2])
     sol = solve(prob)
 
+    @testset "Guesses of initialization problem copied to algebraic variables" begin
+        prob.f.initialization_data.initializeprob[λ] = 1.0
+        prob2 = DiffEqBase.get_updated_symbolic_problem(
+            pend, prob; u0 = prob.u0, p = prob.p)
+        @test prob2[λ] ≈ 1.0
+    end
+
     @testset "`setsym_oop`" begin
         setter = setsym_oop(prob, [Initial(x)])
         (u0, p) = setter(prob, [0.8])

From 997da8d3319ca6b1b238d3988c0ac107fd810dfa Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 28 Apr 2025 11:08:33 +0530
Subject: [PATCH 09/10] fix: handle underdetermined systems in `GetUpdatedU0`

---
 src/systems/problem_utils.jl | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl
index b46ef50f8e..56b7c25f17 100644
--- a/src/systems/problem_utils.jl
+++ b/src/systems/problem_utils.jl
@@ -834,10 +834,13 @@ struct GetUpdatedU0{GA, GIU}
 end
 
 function GetUpdatedU0(sys::AbstractSystem, initsys::AbstractSystem)
-    algevaridxs = BitVector(is_alg_equation.(equations(sys)))
-    algevars = unknowns(sys)[algevaridxs]
+    dvs = unknowns(sys)
+    eqs = equations(sys)
+    algevaridxs = BitVector(is_alg_equation.(eqs))
+    append!(algevaridxs, falses(length(dvs) - length(eqs)))
+    algevars = dvs[algevaridxs]
     get_algevars = getu(initsys, algevars)
-    get_initial_unknowns = getu(sys, Initial.(unknowns(sys)))
+    get_initial_unknowns = getu(sys, Initial.(dvs))
     return GetUpdatedU0(algevaridxs, get_algevars, get_initial_unknowns)
 end
 

From 998a516741251c54fac82114ae25d16aa2a2037c Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 29 Apr 2025 11:36:33 +0530
Subject: [PATCH 10/10] fix: retain explicit initial value, copy guesses for
 unknowns in `GetUpdatedU0`

---
 src/systems/problem_utils.jl | 35 ++++++++++++++++++-----------------
 test/initializationsystem.jl | 11 +++++++++++
 2 files changed, 29 insertions(+), 17 deletions(-)

diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl
index 56b7c25f17..47bf3c678d 100644
--- a/src/systems/problem_utils.jl
+++ b/src/systems/problem_utils.jl
@@ -810,44 +810,45 @@ end
     $(TYPEDEF)
 
 A callable struct to use as the `get_updated_u0` field of `InitializationMetadata`.
-Returns the value of `Initial.(unknowns(sys))`, except with algebraic variables replaced
-by their guess values in the initialization problem.
+Returns the value to use for the `u0` of the problem. 
 
 # Fields
 
 $(TYPEDFIELDS)
 """
-struct GetUpdatedU0{GA, GIU}
+struct GetUpdatedU0{GG, GIU}
     """
-    Mask with length `length(unknowns(sys))` denoting indices of algebraic variables.
+    Mask with length `length(unknowns(sys))` denoting indices of variables which should
+    take the guess value from `initializeprob`.
     """
-    algevars::BitVector
+    guessvars::BitVector
     """
-    Function which returns the values of algebraic variables in `initializeprob`, in the
-    order the algebraic variables occur in `unknowns(sys)`.
+    Function which returns the values of variables in `initializeprob` for which
+    `guessvars` is `true`, in the order they occur in `unknowns(sys)`.
     """
-    get_algevars::GA
+    get_guessvars::GG
     """
     Function which returns `Initial.(unknowns(sys))` as a `Vector`.
     """
     get_initial_unknowns::GIU
 end
 
-function GetUpdatedU0(sys::AbstractSystem, initsys::AbstractSystem)
+function GetUpdatedU0(sys::AbstractSystem, initsys::AbstractSystem, op::AbstractDict)
     dvs = unknowns(sys)
     eqs = equations(sys)
-    algevaridxs = BitVector(is_alg_equation.(eqs))
-    append!(algevaridxs, falses(length(dvs) - length(eqs)))
-    algevars = dvs[algevaridxs]
-    get_algevars = getu(initsys, algevars)
+    guessvars = trues(length(dvs))
+    for (i, var) in enumerate(dvs)
+        guessvars[i] = !isequal(get(op, var, nothing), Initial(var))
+    end
+    get_guessvars = getu(initsys, dvs[guessvars])
     get_initial_unknowns = getu(sys, Initial.(dvs))
-    return GetUpdatedU0(algevaridxs, get_algevars, get_initial_unknowns)
+    return GetUpdatedU0(guessvars, get_guessvars, get_initial_unknowns)
 end
 
 function (guu::GetUpdatedU0)(prob, initprob)
     buffer = guu.get_initial_unknowns(prob)
-    algebuf = view(buffer, guu.algevars)
-    copyto!(algebuf, guu.get_algevars(initprob))
+    algebuf = view(buffer, guu.guessvars)
+    copyto!(algebuf, guu.get_guessvars(initprob))
     return buffer
 end
 
@@ -890,7 +891,7 @@ function maybe_build_initialization_problem(
     initializeprob = remake(initializeprob; p = initp)
 
     get_initial_unknowns = if is_time_dependent(sys)
-        GetUpdatedU0(sys, initializeprob.f.sys)
+        GetUpdatedU0(sys, initializeprob.f.sys, op)
     else
         nothing
     end
diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 2534fb163f..5c36fcba3e 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -1533,6 +1533,17 @@ end
         @test prob2[λ] ≈ 1.0
     end
 
+    @testset "Initial values for algebraic variables are retained" begin
+        prob2 = ODEProblem(
+            pend, [x => (√2 / 2), D(y) => 0.0], (0.0, 1.5),
+            [g => 1], guesses = [λ => 1, y => √2 / 2])
+        sol = solve(prob)
+        @test SciMLBase.successful_retcode(sol)
+        prob3 = DiffEqBase.get_updated_symbolic_problem(
+            pend, prob2; u0 = prob2.u0, p = prob2.p)
+        @test prob3[D(y)] ≈ 0.0
+    end
+
     @testset "`setsym_oop`" begin
         setter = setsym_oop(prob, [Initial(x)])
         (u0, p) = setter(prob, [0.8])