This notebook uses the three point formulas for the approximation to the derivative.

It takes in a function f and approximates it on the interval [a, b] as a column vector of length n+1.
The derivative is approximated as a matrix times the column vector.

a = 0; (* left endpoint *)
b = 2; (* right endpoint *)
n = 5;  (* the intervals [a, b] is broken into n sub-intervals.  We need n ≥ 2 *)
f[x_] :=Exp[2x];
(* Change the stuff above this line *)

h = (b-a)/n;
fList = Table[f[x], {x, a, b, h}];
fpList = Table[f'[x], {x, a, b, h}];
(* The next line defines an n+1 × n+1 sparse matrix (array).  The first pattern that fits is used; for example the (1,2) component of the matrix is 4, event though the Band[] command makes the (i,i+1) component equal to 1 *)
der = 1/(2h) SparseArray[{ {1,1} -> -3, {1,2} -> 4, {1,3} -> -1, {n+1,n+1} -> 3,{n+1, n} -> -4, {n+1,n-1} -> 1,Band[{1,2}]-> 1, Band[{2,1}]-> -1}, {n+1, n+1}];
If[n < 10, Print["der = ", 1/(2h), MatrixForm[2h der] ] ]
ListPlot[{fpList, der.fList}, Joined -> True, DataRange ->{a,b}, PlotLabel ->"Blue is f', orange is approximation to f'"]
ListPlot[ der.fList - fpList, PlotRange ->  All, DataRange ->{a,b}, PlotLabel ->"Approximation to f' - true f'"]