The paper considers numerical differentiation of functions with large gradients in the region of an exponential boundary layer. This topic is important, since the application of classical polynomial difference formulas for derivatives to such functions in the case of a uniform mesh leads to unacceptable errors if the perturbation parameter is comparable with the mesh size. The numerical differentiation formula with a given number of nodes in the difference stencil is built on subintervals covering the original interval. The accuracy of numerical differentiation formulas on a Bakhvalov mesh, which is widely used in the construction of difference schemes for singularly perturbed problems, is analyzed. For the original function of one variable, a representation in the form of a sum of regular and boundary-layer components, based on the Shishkin decomposition, is used to solve a singularly perturbed problem. Previously, such a decomposition was used to prove the convergence of difference schemes. An estimate of the error of classical polynomial formulas for numerical differentiation on a Bakhvalov mesh is obtained. The error estimate on a Bakhvalov mesh is obtained in the general case, when a derivative of an arbitrarily given order is calculated and the difference stencil for this derivative contains a given number of nodes. The error estimate depends on the order of the calculated derivative and the number of nodes in difference stencil and takes into account the uniformity in the parameter . The results of numerical experiments are presented, which are consistent with the error estimates obtained.