`finite_difference_hessian` is mutating arrays

julia> function f(x,p)
       hess = FiniteDiff.finite_difference_hessian(y -> sum(y.^3), x)
       return hess * p
       end
f (generic function with 1 method)

julia> x=rand(3)
3-element Vector{Float64}:
 0.45298015977579165
 0.6696731824795704
 0.8825460798816239

julia> p=rand(3)
3-element Vector{Float64}:
 0.8610782341642329
 0.35801727566906194
 0.8214679367160949

julia> f(x,p)
3-element Vector{Float64}:
 2.3403081413900946
 1.438527409400547
 4.34989985490966

julia> Zygote.gradient(p->sum(f(x,p)), p)[1]
ERROR: Mutating arrays is not supported -- called setindex!(Matrix{Float64}, ...)
This error occurs when you ask Zygote to differentiate operations that change
the elements of arrays in place (e.g. setting values with x .= ...)

Possible fixes:
- avoid mutating operations (preferred)
- or read the documentation and solutions for this error
  https://fluxml.ai/Zygote.jl/dev/limitations.html#Array-mutation-1

Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:33
  [2] _throw_mutation_error(f::Function, args::Matrix{Float64})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\lib\array.jl:70
  [3] (::Zygote.var"#444#445"{Matrix{Float64}})(#unused#::Nothing)
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\lib\array.jl:82
  [4] (::Zygote.var"#2496#back#446"{Zygote.var"#444#445"{Matrix{Float64}}})(Δ::Nothing)
    @ Zygote C:\Users\Luffy\.julia\packages\ZygoteRules\AIbCs\src\adjoint.jl:67
  [5] Pullback
    @ C:\Users\Luffy\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\matmul.jl:514 [inlined]
  [6] Pullback
    @ C:\Users\Luffy\AppData\Local\Programs\Julia-1.7.2\share\julia\stdlib\v1.7\LinearAlgebra\src\matmul.jl:510 [inlined]
  [7] (::typeof(∂(copytri!)))(Δ::Nothing)
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
  [8] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:139 [inlined]
  [9] (::typeof(∂(#finite_difference_hessian!#43)))(Δ::Nothing)
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [10] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:61 [inlined]
 [11] (::typeof(∂(finite_difference_hessian!##kw)))(Δ::Nothing)
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [12] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:41 [inlined]
 [13] (::typeof(∂(#finite_difference_hessian#41)))(Δ::LinearAlgebra.Symmetric{Float64, Matrix{Float64}})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [14] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:39 [inlined]
 [15] (::typeof(∂(finite_difference_hessian##kw)))(Δ::LinearAlgebra.Symmetric{Float64, Matrix{Float64}})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [16] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:31 [inlined]
 [17] (::typeof(∂(#finite_difference_hessian#40)))(Δ::LinearAlgebra.Symmetric{Float64, Matrix{Float64}})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [18] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:30 [inlined]
 [19] (::typeof(∂(finite_difference_hessian)))(Δ::LinearAlgebra.Symmetric{Float64, Matrix{Float64}})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [20] Pullback
    @ C:\Users\Luffy\.julia\packages\FiniteDiff\KkXlb\src\hessians.jl:30 [inlined]
 [21] (::typeof(∂(finite_difference_hessian)))(Δ::LinearAlgebra.Symmetric{Float64, Matrix{Float64}})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [22] Pullback
    @ .\REPL[19]:2 [inlined]
 [23] (::typeof(∂(f)))(Δ::FillArrays.Fill{Float64, 1, Tuple{Base.OneTo{Int64}}})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [24] Pullback
    @ .\REPL[23]:1 [inlined]
 [25] (::typeof(∂(#15)))(Δ::Float64)
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface2.jl:0
 [26] (::Zygote.var"#60#61"{typeof(∂(#15))})(Δ::Float64)
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface.jl:41
 [27] gradient(f::Function, args::Vector{Float64})
    @ Zygote C:\Users\Luffy\.julia\packages\Zygote\IoW2g\src\compiler\interface.jl:76
 [28] top-level scope
    @ REPL[23]:1

Start with src/hessians.jl at the locations shown in the stack trace, especially lines 39, 61, and 139, and reproduce the failure using the Zygote.gradient example from the issue. Done means finite_difference_hessian can be differentiated in this example without triggering Zygote's array-mutation error.