diff --git a/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl b/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl
index 7cae109703..e27001d5d5 100644
--- a/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl
+++ b/lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl
@@ -61,6 +61,33 @@
     return maybe_codegen_scimlfn(expression, NonlinearFunction{iip, spec}, args; kwargs...)
 end
 
+function _get_nonlinear_bounds_arrays(unknowns_list, u0)
+    if !any(hasbounds, unknowns_list)
+        return nothing, nothing
+    end
+
+    T = u0 === nothing ? Float64 : eltype(u0)
+    if !(T <: Number)
+        T = Float64
+    end
+
+    N = length(unknowns_list)
+    lb = Vector{T}(undef, N)
+    ub = Vector{T}(undef, N)
+
+    for (i, v) in enumerate(unknowns_list)
+        b = getbounds(v)
+        if b !== nothing && length(b) >= 2
+            lb[i] = b[1] !== nothing ? T(b[1]) : typemin(T)
+            ub[i] = b[2] !== nothing ? T(b[2]) : typemax(T)
+        else
+            lb[i] = typemin(T)
+            ub[i] = typemax(T)
+        end
+    end
+    return lb, ub
+end
+
 @fallback_iip_specialize function SciMLBase.NonlinearProblem{iip, spec}(
         sys::System, op; expression = Val{false},
         check_length = true, check_compatibility = true, kwargs...
@@ -80,9 +107,24 @@ end
 
     kwargs = process_kwargs(sys; kwargs...)
     ptype = getmetadata(sys, ProblemTypeCtx, StandardNonlinearProblem())
-    args = (; f, u0, p, ptype)
 
-    return maybe_codegen_scimlproblem(expression, NonlinearProblem{_iip}, args; kwargs...)
+    if !haskey(kwargs, :lb) && !haskey(kwargs, :ub)
+        lb, ub = _get_nonlinear_bounds_arrays(unknowns(sys), u0)
+        if lb !== nothing && ub !== nothing
+            kwargs = merge(kwargs, (lb = lb, ub = ub))
+        end
+    end
+
+    if expression == Val{false}
+        prob = SciMLBase.NonlinearProblem(f, u0, p; kwargs...)
+        if !(ptype isa StandardNonlinearProblem)
+            prob = SciMLBase.remake(prob; problem_type = ptype)
+        end
+        return prob
+    else
+        args = ptype isa StandardNonlinearProblem ? (; f, u0, p) : (; f, u0, p, ptype)
+        return maybe_codegen_scimlproblem(expression, NonlinearProblem{_iip}, args; kwargs...)
+    end
 end
 
 @fallback_iip_specialize function SciMLBase.NonlinearLeastSquaresProblem{iip, spec}(
@@ -100,11 +142,22 @@ end
     )
 
     kwargs = process_kwargs(sys; kwargs...)
-    args = (; f, u0, p)
 
-    return maybe_codegen_scimlproblem(
-        expression, NonlinearLeastSquaresProblem{_iip}, args; kwargs...
-    )
+    if !haskey(kwargs, :lb) && !haskey(kwargs, :ub)
+        lb, ub = _get_nonlinear_bounds_arrays(unknowns(sys), u0)
+        if lb !== nothing && ub !== nothing
+            kwargs = merge(kwargs, (lb = lb, ub = ub))
+        end
+    end
+
+    if expression == Val{false}
+        return SciMLBase.NonlinearLeastSquaresProblem(f, u0, p; kwargs...)
+    else
+        args = (; f, u0, p)
+        return maybe_codegen_scimlproblem(
+            expression, NonlinearLeastSquaresProblem{_iip}, args; kwargs...
+        )
+    end
 end
 
 function check_compatible_system(
diff --git a/lib/ModelingToolkitBase/test/nonlinearsystem.jl b/lib/ModelingToolkitBase/test/nonlinearsystem.jl
index b97d352774..462355c70c 100644
--- a/lib/ModelingToolkitBase/test/nonlinearsystem.jl
+++ b/lib/ModelingToolkitBase/test/nonlinearsystem.jl
@@ -515,3 +515,42 @@ end
     prob = NonlinearProblem(sys, [x => 1.0])
     @test prob.problem_type == "A"
 end
+
+@testset "Bounds extraction for Nonlinear Problems" begin
+    @variables x [bounds=(0.0, 1.0)]
+    @variables y [bounds=(-1.0, 5.0)]
+    @variables z
+
+    eqs = [
+        0 ~ x^2 - 0.5,
+        0 ~ y - x,
+        0 ~ z - y
+    ]
+
+    @named sys = System(eqs, [x, y, z], [])
+    sys = complete(sys)
+
+    u0 = [x => 0.5, y => 0.5, z => 0.5]
+
+    prob_nl = NonlinearProblem(sys, u0)
+    prob_ls = NonlinearLeastSquaresProblem(sys, u0)
+
+    for p in (prob_nl, prob_ls)
+        @test p.lb !== nothing
+        @test p.ub !== nothing
+
+        @test eltype(p.lb) === Float64
+        @test eltype(p.ub) === Float64
+
+        unks = unknowns(sys)
+        expected_lb = map(unks) do v
+            isequal(v, x) ? 0.0 : (isequal(v, y) ? -1.0 : typemin(Float64))
+        end
+        expected_ub = map(unks) do v
+            isequal(v, x) ? 1.0 : (isequal(v, y) ? 5.0 : typemax(Float64))
+        end
+
+        @test p.lb == expected_lb
+        @test p.ub == expected_ub
+    end
+end