function df=dnf(f,x,n,dx)
% degree nth derivative of a function f
if nargin<4  
    dx=0.001;
end
if nargin<3 
    dx=0.001;n=1;
end
dxh=1.0/dx;
if n==0 df=f(x);
elseif n==1
df=(3.0*f(x-4.0*dx)-32.0*f(x-3.0*dx)+168.0*f(x-2.0*dx)-672.0*f(x-dx)+672.0*f(x+dx)-168.0*f(x+2.0*dx)+32.0*f(x+3.0*dx)-3.0*f(x+4.0*dx))/840.0*dxh;
elseif n==2
df=(-14350.0*f(x)-9.0*f(x-4.0*dx)+128*f(x-3.0*dx)-1008*f(x-2.0*dx)+8064*f(x-dx)+8064.0*f(x+dx)-1008.0*f(x+2.0*dx)+128.0*f(x+3.0*dx)-9.0*f(x+4.0*dx))/5040.0*dxh*dxh;
elseif n==3
df=(-7.0*f(x-4.0*dx)+72.0*f(x-3.0*dx)-338.0*f(x-2.0*dx)+488.0*f(x-dx)-488.0*f(x+dx)+338.0*f(x+2.0*dx)-72.0*f(x+3.0*dx)+7.0*f(x+4.0*dx))/240.0*dxh*dxh*dxh;
elseif n==4
df=(2730.0*f(x)+7.0*f(x-4.0*dx)-96.0*f(x-3.0*dx)+676.0*f(x-2*dx)-1952.0*f(x-dx)-1952.0*f(x+dx)+676.0*f(x+2.0*dx)-96.0*f(x+3.0*dx)+7.0*f(x+4.0*dx))/240.0*dxh*dxh*dxh*dxh;
elseif n==5
df=(f(x-4.0*dx)-9.0*f(x-3.0*dx)+26.0*f(x-2.0*dx)-29.0*f(x-dx)+29.0*f(x+dx)-26.0*f(x+2.0*dx)+9.0*f(x+3.0*dx)-f(x+4.0*dx))/6.0*dxh*dxh*dxh*dxh*dxh;
elseif n==6
df=(-150.0*f(x)-f(x-4.0*dx)+12.0*f(x-3.0*dx)-52.0*f(x-2.0*dx)+116.0*f(x-dx)+116.0*f(x+dx)-52.0*f(x+2.0*dx)+12.0*f(x+3.0*dx)-f(x+4.0*dx))/4.0*dxh*dxh*dxh*dxh*dxh*dxh;
elseif n==7
df=(-f(x-4.0*dx)+6.0*f(x-3.0*dx)-14.0*f(x-2.0*dx)+14.0*f(x-dx)-14.0*f(x+dx)+14.0*f(x+2.0*dx)-6.0*f(x+3.0*dx)+f(x+4.0*dx))/2.0*dxh*dxh*dxh*dxh*dxh*dxh*dxh;
elseif n==8
df=(70.0*f(x)+f(x-4.0*dx)-8.0*f(x-3.0*dx)+28.0*f(x-2.0*dx)-56.0*f(x-dx)-56.0*f(x+dx)+28.0*f(x+2.0*dx)-8.0*f(x+3.0*dx)+f(x+4.0*dx))*dxh*dxh*dxh*dxh*dxh*dxh*dxh*dxh;

else df=0;end
end