Guard JacobianCache callable for empty u, add empty u0 tests

Also guard the JacobianCache callable to early-return the empty Jacobian
when length(u) == 0, preventing DI.jacobian! from being called with
nothing di_extras.

Add comprehensive test for the u0=nothing pathway: IIP and OOP
NonlinearProblems with Float64[] u0, SymbolCache-based sys with
observed quantities, parameter mutation via symbolic indexing, direct
solve, caching interface, and step interface.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: ChrisRackauckas

lib/NonlinearSolveBase/src/jacobian.jl

@@ -176,6 +176,9 @@ end
 
 ## Actually Compute the Jacobian
 function (cache::JacobianCache)(u)
+    # Empty state vector — Jacobian is 0×0, nothing to compute
+    length(u) == 0 && return cache.J
+
     cache.stats.njacs += 1
     (; f, J, p) = cache
     if SciMLBase.isinplace(f)

test/core_tests.jl

@@ -542,3 +542,71 @@ end
     solve!(iter)
     @test precs.i == 1
 end
+
+@testitem "Empty u0 (u0=nothing pathway)" tags = [:core] begin
+    using SymbolicIndexingInterface
+
+    # Mimics MTK systems with no state variables — only parameters and observed quantities.
+    # OrdinaryDiffEqCore converts u0=nothing to Float64[] before reaching NonlinearSolve.
+    sys = SymbolCache(nothing, Dict(:a => 1, :b => 2), nothing)
+
+    function observed_fn(sym)
+        if sym == :result
+            return (u, p) -> p[1] + p[2]
+        end
+        return nothing
+    end
+
+    @testset "IIP" begin
+        f_iip!(resid, u, p) = nothing
+        nlfunc = NonlinearFunction(f_iip!; observed = observed_fn, sys = sys)
+        prob = NonlinearProblem{true}(nlfunc, Float64[], [3.0, 4.0])
+
+        # init/solve cycle
+        cache = @test_nowarn init(prob, NewtonRaphson())
+        @test state_values(cache) == Float64[]
+        @test parameter_values(cache) == [3.0, 4.0]
+
+        # parameter mutation through symbolic indexing
+        cache.ps[:a] = 10.0
+        cache.ps[:b] = 20.0
+        @test parameter_values(cache) == [10.0, 20.0]
+
+        sol = solve!(cache)
+        @test SciMLBase.successful_retcode(sol)
+        @test isempty(sol.u)
+    end
+
+    @testset "OOP" begin
+        f_oop(u, p) = Float64[]
+        nlfunc = NonlinearFunction(f_oop; observed = observed_fn, sys = sys)
+        prob = NonlinearProblem(nlfunc, Float64[], [5.0, 6.0])
+
+        sol = @test_nowarn solve(prob, NewtonRaphson())
+        @test SciMLBase.successful_retcode(sol)
+        @test isempty(sol.u)
+    end
+
+    @testset "Default algorithm" begin
+        f_iip!(resid, u, p) = nothing
+        nlfunc = NonlinearFunction(f_iip!; observed = observed_fn, sys = sys)
+        prob = NonlinearProblem{true}(nlfunc, Float64[], [1.0, 2.0])
+
+        sol = @test_nowarn solve(prob)
+        @test SciMLBase.successful_retcode(sol)
+        @test isempty(sol.u)
+    end
+
+    @testset "Step interface" begin
+        f_iip!(resid, u, p) = nothing
+        nlfunc = NonlinearFunction(f_iip!; observed = observed_fn, sys = sys)
+        prob = NonlinearProblem{true}(nlfunc, Float64[], [1.0, 2.0])
+
+        cache = init(prob, NewtonRaphson())
+        for i in 1:10
+            step!(cache)
+            cache.force_stop && break
+        end
+        @test SciMLBase.successful_retcode(cache.retcode)
+    end
+end