Via Taylor's theorem

We combine Taylor approximations of $f$ at different points “in the right way” to cancel error terms and obtain an approximation of the function derivative.
This technique is perhaps best illustrated by example.
If $f$ has two continuous derivatives on $(a,b)$, then Taylor's theorem guarantees that for every $x_0 \in(a,b)$ there exists $x_0 < \xi < x_0 + h$ such that $$f(x_0+ h) = f(x_0) + h \color{var(--emphColor)}{f'(x_0)} +\frac{h^2}{2}f''(\xi).$$
Simply isolating $\color{var(--emphColor)}{f'(x_0)}$ defines an approximation formula: \begin{equation} \color{var(--emphColor)}{f'(x_0)} = \frac{f(x_0 + h) - f(x_0)}{h}- \frac{h}{2} f''(\xi). \end{equation}
We assume $h > 0$, so the formula $$\color{var(--emphColor)}{f'(x_0)} \approx \frac{f(x_0 + h) - f(x_0)}{h}$$ “steps forward” to approximate the derivative, so this formula defines a forward difference approximation.
In addition, the error term $\frac{h}{2} f''(\xi)$ is proportional to $h$, so our approximation is first order.
In practice, when the error is proportional to $h$, we expect that halving the step size approximately cuts the error in half.
What if we want the error to decrease faster?
We'll use a Taylor table to derive higher order formulas!
Set $f_{j+k} = f(x_j + kh).$
Suppose we want to use $f_{j-1}, f_j$, and $f_{j+1}$ to obtain a second order formula approximating $f'(x_j)$.
Then we must find coefficients $\alpha_{-1}, \alpha_0$, and $\alpha_1$ such that $$f'(x_j) = \alpha_{-1}f_{j-1} + \alpha_0 f_j +\alpha_{-1} f_{j+1} + O(h^\color{var(--emphColor)}{2}).$$
Taylor's theorem implies that \begin{equation} f_{j+k} = f(x_j) + khf'(x_j)+ \frac{(kh)^2}{2} f''(x_j) + O(h^3). \end{equation}
Substituting into our desired expression, we obtain \begin{align} f'(x_j)& = \alpha_{-1} f_{j\color{var(--emphColor)}{-}1} + \alpha_0 f_j + \alpha_1 f_{j\color{var(--emphColor)}{+}1} + O(h^2)\\ & = \alpha_{-1}\big(f(x_j) \color{var(--emphColor)}{-} h f'(x_j) +\frac{h^2}{2} f''(x_j) + O(h^3)\big)\\ & + \alpha_0 f(x_j)\\ & + \alpha_1 \big(f(x_j) \color{var(--emphColor)}{+} hf'(x_j) + \frac{h^2}{2} f'' (x_j) + O(h^3)\big) \end{align}
Now we combine like terms and write \begin{align} f'(x_j) &= (\alpha_{-1} + \alpha_0 + \alpha_1)f(x_j)\\ & \,\,\, + (\color{var(--emphColor)}{-}h\alpha_{-1} \color{var(--emphColor)}{+} h \alpha_1) f'(x_j)\\ & \,\,\, + \big(\frac{h^2}{2} \alpha_{-1} + \frac{h^2}{2} \alpha_1\big) f'' (x_j)\\ & \,\,\, + (\alpha_{-1} + \alpha_1) O(h^3). \end{align}
Finally, we compare coefficients on both sides to obtain a system of equations.
  • The coefficient of $f(x_j)$ on the LHS is $0$. On the RHS, it is $\alpha_{-1} + \alpha_0 + \alpha_1$, so we must have $$\alpha_{-1} + \alpha_0 + \alpha_1 = 0.$$
  • The coefficient of $f'(x_j)$ on the LHS is $1$. On the RHS, it is $-h\alpha_{-1} + h \alpha_1$, so we must have $$-h\alpha_{-1} + h \alpha_1 = 1.$$
  • The term $f''(x_j)$ does not appear on the LHS, and its coefficient on the RHS is $\frac{h^2}{2} \alpha_{-1} + \frac{h^2}{2} \alpha_1$, so we must have $$\frac{h^2}{2} \alpha_{-1} + \frac{h^2}{2} \alpha_1 = 0.$$
(Continued.) As a single matrix equation, our system reads \begin{align} \begin{bmatrix} 1 & 1 & 1\\ -1 & 0 & 1\\ 1/2 & 0 & 1/2 \end{bmatrix} \begin{bmatrix} \alpha_{-1}\\ \alpha_0\\ \alpha_1 \end{bmatrix}= \frac{1}{h} \begin{bmatrix} 0\\ 1\\ 0 \end{bmatrix} \end{align}
We can input $A$ and $e_2$ into MATLAB to find $x= A\backslash e_2$. We may then set $$\begin{bmatrix} \alpha_{-1}\\ \alpha_0\\ \alpha_1 \end{bmatrix} = \frac{1}{h}x.$$
Here $x$ has numerical values and $h$ is a symbol we use to write our formula.
(Continued.) As a single matrix equation, our system reads \begin{align} \begin{bmatrix} 1 & 1 & 1\\ -1 & 0 & 1\\ 1/2 & 0 & 1/2 \end{bmatrix} \begin{bmatrix} \alpha_{-1}\\ \alpha_0\\ \alpha_1 \end{bmatrix}= \frac{1}{h} \begin{bmatrix} 0\\ 1\\ 0 \end{bmatrix}. \end{align} In this case we don't need MATLAB!