function showgrad(fun,x0,par)

% Gradient Checker
%
% fun = function name
% x0  = point the check is done in (vector or matrix)
% par = parameters passed to the function
% 
% [f,dfdx] = fun(x,par)

vecsize = length(x0);
[xrow,xcol] = size(x0);
funcall = sprintf('[f,g]=%s(x,par);',fun);
xvec=linspace(-0.1,0.1,51); 
dx=xvec(2)-xvec(1);

for ii=1:xrow
    for jj=1:xcol
        
        dataf=[];
        datag=[];
        
        for del_x=xvec,
            x = x0;
            x(ii,jj) = x0(ii,jj) + del_x;
            eval(funcall);
            dataf = [dataf;f];
            datag = [datag;g(ii,jj)];
        end
        
        difff = diff(dataf);
        datagnum = [difff(1);difff]/dx;
        
        plot(xvec,datag,'-',xvec,datagnum,'.')
        title(sprintf('Coordinate i=%d, j=%d',ii,jj))
        legend('Analytical','Numerical')
        if ~((ii+jj)==(xrow+xcol)), pause; end
        
    end
end