Calculating Third (and Higher) Derivatives: ReverseDiff gives 0s while ForwardDiff works

I am trying to calculate the third derivative (tensor). ForwardDiff seems to work but is incredibly slow, also for lower orders. ReverseDiff, while fast for lower derivatives, doesn’t seem to work for third derivatives or higher. Below is a simple example.

using ForwardDiff
using ReverseDiff

x = [1. 3.]
f = z -> cos(z[1])*cos(z[2])
println(ForwardDiff.jacobian(xh -> ForwardDiff.hessian(f, xh), x))
println(ReverseDiff.jacobian(xh -> ReverseDiff.hessian(f, xh), x))
[-0.833049961066805 0.07624746575887673; 0.07624746575887673 -0.833049961066805; 0.07624746575887673 -0.833049961066805; -0.833049961066805 0.07624746575887673]
[0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]

I will try TaylorSeries next, as suggested in another post for higher order derivative computation. In the meantime, I wanted to confirm wether I am using ReverseDiff incorrectly or if there is another issue.

I would try TaylorDiff, which is designed for this kind of problem.

I would also suggest x = [1, 3] rather than [1 3] (that is, use a 1d array, a vector, not a row vector).