I am implementing on Matlab a high-order finite differences scheme to approximate the first derivative of $f(x_i)$ given $x = [x(1), x(2),..., x(i),..., x(n)]$ and $f = [f(x(1)),..,f(x(n))]$ with $x$ a non-uniform grid from $x(1)$ to $x(n)$ (in fact they are Legendre-Gauss-Lobatto nodes in $[-1,1]$). Right now I use the method detailed in chapter 2.1 at https://en.wikibooks.org/wiki/Using_High_Order_Finite_Differences/Definitions_and_Basics#approximation_of_a_first_derivative, which uses the Vandermonde matrix.
I decided to first set $h = x(i+1)-x(i)$ or $x(i)-x(i-1)$ depending on $i$, then I compute the Vandermonde matrix $M$ using $\alpha(j) = (x(j)-x(i))/h$ and I get the $a_j$'s coefficients with $M\backslash [0,1,0,..,0]$'. And finally $f'(x_i) ~ (1/h)*(a_1*f(x_1) + ... + a_n*f(x_n))$. It seems like this method returns the right result, unfortunately it makes the computations really slow...
Do you know another method to compute arbitrary order finite differences schemes ? Or a faster way to compute the $a_j$'s coefficients ?
I hope my explanations are clear enough, and already thank you for your answers.