## Author: M. Turhan ÇOBAN
## Created: 2021-06-06

function retval = dfunc (f,x,x_ref)
      #derivative of the function with respect to x_ref
      h0=0.015675863;
      m=7;
      n=length(x);
      for i=1:n
          x1(i)=x(i);
          x2(i)=x(i);
      end
      #derivative of a simple function
      T=zeros(m,m);
      for i=1:m
          h(i)=0;
          for j=1:m
              T(i,j)=0;
          end
      end
      h(1)=h0;
      r=0.5;
      for i=2:m
          h(i)=h0*r^i;
      end
      for i=1:m
          x1(x_ref)=x1(x_ref)+h(i);
          x2(x_ref)=x2(x_ref)-h(i);
          f1=f(x1);
          f2=f(x2);
          T(i,1)=( f1 - f2)/(2.0*h(i));
          x1(x_ref)=x(x_ref);
          x2(x_ref)=x(x_ref);
      end
      for k=2:m
           for i=1:m-k
               T(i,k)=(h(i)*h(i)*T(i+1,k-1) - h(i+k)*h(i+k)*T(i,k-1))/(h(i)*h(i)- h(i+k)*h(i+k));
           end
      end
      retval=T(1,n-1);
endfunction
