In many nonlinear problems it is necessary to estimate the Jacobian matrix of a nonlinear mapping $F$. In large scale problems the Jacobian of $F$ is usually sparse, and then estimation by differences is attractive because the number of differences can be small compared to the dimension of the problem. For example, if the Jacobian matrix is banded then the number of differences needed to estimate the Jacobian matrix is, at most, the width of the band. In this paper we describe a set of subroutines whose purpose is to estimate the Jacobian matrix of a mapping $F$ with the least possible number of function evaluations.