function res = gradfun(fun,x0,varargin)
% gradfun(fun,x,params)
% gradient of a function at position x0
seps = sqrt(eps); %! Make the step small
h = x0.*seps; %! in each direction 
h(abs(x0)<seps) = seps; %! But not too small!
f0 = feval(fun,x0,varargin{:}); %! Find $f({\bf x}_0)$.
res = zeros(size(x0));
for k=1:length(x0) %! Consider each direction
    x1 = x0;
    x1(k) = x0(k)+h(k); %! Displaced in that direction
    f1 = feval(fun,x1,varargin{:});
    res(k) = (f1 - f0)./h(k); %! The partial derivative
end