Hessian-vector product with central finite differences

[bob-carpenter]

Description

We are implementing a Hessian-vector product for BridgeStan (roualdes/bridgestan#126) and would like to be able to evaluate with central finite differences for cases where forward-mode is not defined (mainly implicit functions like the ode solvers).

There is also a scheme for setting relative perturbations that @justindomke mentioned in the issue above that we should try. Code sketches and reference links are in the above-referenced BridgeStan issue.

Expected Output

This should be a functor defined parallel to where Hessians are defined via finite differences.

Current Version:

v4.6.2