We'll focus on the following two questions.

Given a smooth function $f$ defined by a data table, how can I approximate the value of $f'$ at a point?
Given a function $f$, or even just an oracle evaluating $f$, how can I approximate values of $f'$ by linear combinations of values of $f$?
This is common in applications dealing with complex functions, like those describing forces in molecular dynamics simulations, or deep neural networks.
To be fair, ML frameworks like TensorFlow also provide oracles for evaluating gradients analytically.

There are two main approaches.

Given a data table, construct the interpolating polynomial (or spline) and evaluate the polynomial's derivative analytically.
The exact derivative of the polynomial interpolant serves to approximate the derivative of the function defined by the data.
Use a linear system!
Use Taylor's theorem to find which linear combinations of function values approximate the function derivative at a given point.

Recall that the Lagrange form of the polynomial interpolating the data points $(x_0, f_0), \ldots,(x_n,f_n)$ is given by $$p(x) =\sum^n_{j=0} f_jL_{j}(x).$$

Setting $d_{ij} = L'_{j}(x_i)$, we have \begin{align} \begin{bmatrix} f'(x_0)\\ \vdots\\ f'(x_n) \end{bmatrix} \color{var(--emphColor)}{\approx} \begin{bmatrix} p'(x_0)\\ \vdots\\ p'(x_n) \end{bmatrix} = \begin{bmatrix} \sum_j d_{0,j}f_j\\ \vdots\\ \sum_j d_{n,j}f_j \end{bmatrix} =D\mathbf{f}, \end{align} with $\mathbf{f} = [f_0\cdots f_n]^T$.

A direct calculation shows that \begin{align} d_{ij} = L'_{n,j}(x_i)= \begin{cases} \frac{w_j/w_i}{x_i - x_j}, & i\not=j\\ \sum_{i\not=j} L'_n,j(x_i), & i =j. \end{cases} \end{align}