function dertable=ndre(Fun,x0,h,M,n);
%%  function dertable=ndre(Fun,x0,h,M,n);
%%
%%  Numerical differentiation using the formula:
%%       f'(x_0)=(f(x0+h)-f(x0-h))/2h
%%  and Richardson Extrapolation
%%
%%  INPUT:
%%    Fun - the function data which may be entered EITHER as
%%           a string, enclosed in ' ', with the name the m-file of the 
%%           function ,  
%%           OR as a row vector of data [N1(h),N1(h/2),...,N1(h/2^M)]
%%     x0 -  the point at which the derivative is to be estimated
%%     h  -  the step size h (several h may be given in a row vector
%%            when Fun is given by a m-file).
%%     M  -  maximum ``depth'', that is, h/2^M
%%     n  -  maximum movement sideways (number of extrapolations)
%%  OUTPUT:
%%    dertable - table of approximate values of the derivative 


if isstr(Fun)==1;  % test whether from a function file or from given data
   dertable(1,1)=nd4_13(Fun,x0,h);
   for j=1:M;
     h=h/2;
     dertable(j+1,1)=nd4_13(Fun,x0,h);
     for k=1:min(j,n);
        dertable(j+1,k+1)=(4^k*dertable(j+1,k)-dertable(j,k))/(4^k-1);
     end;
   end;
else
   dertable(:,1)=Fun';
   for j=1:M;
     for k=1:max(j,n);
        dertable(j+1,k+1)=(4^k*dertable(j+1,k)-dertable(j,k))/(4^k-1);
     end;
   end;
end;
