Given a function $f: [x_0,x_4] → \Bbb R$ and equidistant points $x_0, x_1, x_2, x_3, x_4$ so that $h=x_{i+1} - x_i > 0$.
Normally I would do, $f'(x) \approx \dfrac{f(x+h)-f(x)}{h}$, but here I want to find approximation for $f'(x_2)$ as linear combination of function values in these points as $f'(x_2) \approx \sum_{j=0}^4 \alpha_j f(x_j)$. There are no boundary conditions so I can't make a system of linear equations, at least I see no straight forward way.
Are there any standard methods of how to determine alphas related to the finite elements method, so that the formula will be as precise as possible for polynomials of as high order as possible?
I see that this is closely related to the five-point midpoint formula https://www3.nd.edu/~zxu2/acms40390F15/Lec-4.1.pdf if one sets x_0= x_2 -2h and so on, but how to derive this formula?
How does that generalise to higher order derivatives?
Use Taylor series for small values of $h$ $$f(x)=f(x+(k-2) h)=f(x)+h (k-2) f'(x)+\frac{1}{2} h^2 (k-2)^2 f''(x)+\frac{1}{6} h^3 (k-2)^3 f^{(3)}(x)+\frac{1}{24} h^4 (k-2)^4 f^{(4)}(x)+\frac{1}{120} h^5 (k-2)^5 f^{(5)}(x)+O\left(h^6\right)$$
Now, write $$h f'(x)=a f(x-2h)+b f(x-h)+c f(x)+d f(x+h)+e f(x+2h)$$ and replace to get $$h f'(x)=f(x) (a+b+c+d+e)+h f'(x) (-2 a-b+d+2 e)+\frac{1}{2} h^2 f''(x) (4 a+b+d+4 e)+\frac{1}{6} h^3 f^{(3)}(x) (-8 a-b+d+8 e)+\frac{1}{24} h^4 f^{(4)}(x) (16 a+b+d+16 e)+\frac{1}{120} h^5 f^{(5)}(x) (-32 a-b+d+32 e)+O\left(h^6\right)$$ Then, you have the equations $$a+b+c+d+e=0 \tag 1$$ $$-2 a-b+d+2 e=1 \tag 2$$ $$4 a+b+d+4 e=0 \tag 3$$ $$-8 a-b+d+8 e=0 \tag 4$$ $$16 a+b+d+16 e=0 \tag 5$$
from which $a=\frac 1 {12}$, $b=-\frac 2 {3}$, $c=0$, $d=\frac 2 {3}$, $e=-\frac 1 {12}$.
So, using your notations $$f'(x_2)=\frac{f(x_0)-8f(x_1)+8f(x_3)-f(x_4) }{12 h}$$ and the "error" is $\frac 1 {30} h^4 f^{(5)}(x_2)$
Let us try for $$f(x)=e^x \tan (x) \log (\sin (x)+1)$$ at $x=2$ using $h=0.1$. This would give for the derivative $f'(2)=20.671695$ while the exact value would be $20.671717$