The solution of a nonlinear optimization problem often requires an estimate of the Hessian matrix for a function $f$. In large scale problems the Hessian matrix is usually sparse, and then estimation by differences of gradients is attractive because the number of differences can be small compared to the dimension of the problem. In this paper we describe a set of subroutines whose purpose is to estimate the Hessian matrix with the least possible number of gradient evaluations.