Support immutable output types in out-of-place scalar→vector gradient

[ChrisRackauckas-Claude]

Summary

• finite_difference_gradient (out-of-place, cached) for scalar x previously delegated to finite_difference_gradient!, which uses @. df = result / epsilon. This in-place broadcast calls copyto! → setindex! on the output buffer, which fails for immutable array types like SVector or ArrayPartition{SVector}.
• Extract the scalar→vector case into _scalar_gradient_oop which computes the finite difference purely out-of-place using @. (a - b) / h (allocating a new array via copy rather than mutating via copyto!).
• Add regression test using a ReadOnlyVec wrapper that blocks setindex! to simulate immutable array types.

Motivation

When solving a SecondOrderODEProblem with SVector state in OrdinaryDiffEq.jl, the solver creates an ArrayPartition(SVector, SVector) state. If the solver switches to Rosenbrock23 with finite differences, calc_tderivative calls DI.derivative → finite_difference_gradient, which allocates a buffer via zero(cache.c1). Since similar(::ArrayPartition{SVector}) preserves the immutable SVector inner parts (unlike similar(::SVector) which returns MVector), the in-place @. df = ... fails with setindex!(::SVector{N,T}, ...) is not defined.

The out-of-place finite_difference_gradient API should never require mutation of the result — it already allocates a return value. This PR makes it compute directly without mutation for the scalar input case.

Fixes SciML/OrdinaryDiffEq.jl#3444

Test plan

• Existing tests pass (all 161 passed, 1 pre-existing broken)
• New test verifies finite_difference_gradient works with immutable output arrays (forward and central)
• Verified fix resolves the original OrdinaryDiffEq.jl issue: SecondOrderODEProblem with SVector state + Rosenbrock23(autodiff=false) now solves successfully

🤖 Generated with Claude Code