function p=df(f,x,dx)
% derivative of a function
if nargin<3 
    dx=0.001;
end
p = ( -3*f(x+4*dx)+32*f(x+3*dx)-168*f(x+2*dx)+672*f(x+dx)-672*f(x-dx) + 168*f(x-2*dx)-32*f(x-3*dx)+3*f(x-4*dx) )/(dx*840); 
end