Efficient Arbitrary Order Finite Differences in 1D

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 ?