Numerical Hessian (Multivariate)

["Anatoliy Belaygorod" <belaygorod at wustl dot edu>]

Hello,
Can somebody please tell me how I can use GSL to numerically evaluate a high-dimensional function's Hessian at some point?
 
I know that GSL has "gsl_diff_central " that can be used to evaluate a GRADIENT of a high-dimensional function at some point simply by looping through all the parameters and differentiating one at the time. Can something like this be used to evaluate a Hessian?
 
I suspect that GSL has UNDOCUMENTED Hessian function, because it is needed internally for a variety of routines implemented in GSl, e.g. multidimensional minimization most definitely made use of some numerical Hessian routine. Unfortunately, it seems not to be documented.
 
Any code or example would be HIGHLY appreciated.
 
Sincerely,
Anatoliy