Can a neural network learned function $f$ be used for calculating its derivate $f'$? $y$ is a scalar value, $x$ is an $m \times 1$ vector. Function $f$ is a map $y = f(x)$.
Say we have $n$ examples {$(y_1, x_1), (y_2, x_2), \cdots, (y_n, x_n)$}, and we use these examples for fitting a feed-forward neural network. After the NN is optimized, we get a proxy $f^*$ of the underlying function $f$. 
Can this $f^*$ be used for calculating the derivate of $f$? 
How can we calculate the derivate using this $f^*$? 
Say we do this by following the definition of the derivate of a function 
$$f' = \text{lim}_{\Delta x \rightarrow 0} \frac{f(x + \Delta x ) - f(x)}{\Delta x}$$
However, when $\Delta x$ is a vector, how can we choose this $\Delta x$ for calculating the limit? 
 A: While you can use finite differences, there is a better way: Just as reverse mode automatic differentiation is used for backpropagation, you can also use it for computing the gradient of $f$ wrt the inputs. This is faster and gives an exact answer.
A: As @shimao suggested, ideally you implement the neural network in one of the many symbolic computation frameworks (such as Theano or Tensor Flow). They have the great benefit that they represent the neural network as an analytical formula and can automatically compute its analytical derivatives.
You can also perform some finite differences numerical approximation using the neural network, but it will be computationally more expensive. One thing that you need to take into account that differentiation of a function taking vector input is not a single number but gradient, a vector of derivatives w.r.t. elements of the input vector:
$$
\nabla f=\left(\frac{\partial f}{\partial x_0}, \ldots, \frac{\partial f}{\partial x_m}\right),
$$
where $\frac{\partial f}{\partial x_i}$ is defined as
$$
\frac{\partial f}{\partial x_i} = \text{lim}_{\varepsilon \rightarrow 0} \frac{f(x + \varepsilon \mathbb{I}_i ) - f(x)}{\varepsilon}
$$
and $\mathbb{I}_i$ is an indicator vector whose $i$-th element is one and other elements are zeros. Computing gradient with this method will require $m+1$ function evaluations, as opposed to taking its analytical representation which is a single function evaluation.
Now regarding the actual value of $\varepsilon$ to use, it should be a tiny number. For example, the numerical differentiation implemented in scipy uses square root of machine epsilon for float64 data type, referencing [1] as the source for this choice.
Finally, note that $f^*$ is only an approximation of $f$ and its derivatives do not necessarily have to be good approximations of the derivatives of $f$, especially in the parts of the input space that are not covered well in the training data.

[1] W. H. Press et. al. "Numerical Recipes. The Art of Scientific
           Computing. 3rd edition", sec. 5.7.
