function p=df(f,x,dx)
% derivative of a function 
if nargin == 2 dx=0.001;
p=-3*f(x+4*dx)+32*f(x+3*dx)-168*f(x+    
iter=0;
es=0.000001;
ea=1.1*es;
ia=0;
iu=0;
fa=f(xa);
fu=f(xu);
p=xu;
double xpold;
while(iter<maxit && ea>es)
xpold=p;    
p=xu-fu*(xa-xu)/(fa-fu);
fp=f(p);
iter=iter+1;
if p~=0 ea=abs((p-xpold)/p)*100.0;
test=fa*fp;
if test<0 
    xu=p;
    fu=f(xu);
    iu=0;ia=ia+1;
    if ia >= 2 fa=fa/2.0;end
elseif test>0  xa=p;fa=f(xa);ia=0;iu=iu+1;
        if iu>=2 fu=fu/2.0; end
elseif test==0 ea = 0; end
end
if iter>=maxit fprintf('Maximum number of iteration is exceeded result might not be valid'); end
end
