function res = derivfun(fun,x0,varargin)
% derivfun(fun,x0,params)
% numer. derivative of a function at position x0
seps = sqrt(eps); %! make the step small
h = x0.*seps;
h(abs(x0)<seps) = seps; %! h when x0 is close to zero
x1 = x0 + h;
f0 = feval(fun,x0,varargin{:});
f1 = feval(fun,x1,varargin{:});
res = (f1-f0)./(x1-x0);