function res = jacobian(fun,x0,varargin)
% jacobian(fun,x,params)
% jacobian of a vector valued function at position x0
seps = sqrt(eps); %! Make the step small
h = x0.*seps; %! a vector 
h(abs(x0)<seps) = seps; %! But not too small!
displaced = x0 + h; %! A step in all directions
stepsize = displaced- x0;
f0 = feval(fun,x0,varargin{:});
res = zeros(length(f0(:)),length(x0(:)));
for k=1:length(x0(:)) %! Consider each direction
    x1 = x0;
    x1(k) = displaced(k); %! Displaced in that direction
    f1 = feval(fun,x1,varargin{:});
    partial = (f1 - f0)./stepsize(k); %! The partial derivative
    res(:,k) = partial(:);
end