function [df]=numgrad(func,x0,method,varargin);
%
% Numerical gradient
%
% [df]=numgrad(func,x0,method,param)
%
% method = 'c' -> centered
%        = 'l' -> left
%        = 'r' -> right	
%
x0		= x0(:);
f		= feval(func,x0,varargin{:});
m		= length(x0);
n		= length(f);
df		= zeros(n,m);
dev	= diag(.00001*max(abs(x0),1e-8*ones(size(x0))));
if (lower(method)=='l');
   for i=1:m;
      df(:,i)	= (feval(func,x0+dev(:,i),varargin{:})-f)/dev(i,i);
   end;
elseif (lower(method)=='r')
   for i=1:m;
      df(:,i)	= (f-feval(func,x0-dev(:,i),varargin{:}))/dev(i,i);
   end;
elseif (lower(method)=='c')
   for i=1:m;
      df(:,i)	= (feval(func,x0+dev(:,i),varargin{:})-feval(func,x0-dev(:,i),varargin{:}))/(2*dev(i,i));
   end;
else
   disp('Bad method specified')
end
