function [f,g] = gradx(func,x,der)
%
% gradient computation by 1st/2nd forward or central finite difference for a function
%
if nargin < 3;der = 2;end
if der < 3
    h = sqrt(eps);
else
    h = eps^(1/3);
end
n = length(x);
f = func(x);
g = zeros(n,1);

for i = 1:n
   xp = x;
   xp(i) = xp(i) + h;
   
   switch der
       case 1
       g(i)= (func(xp) -f)/h;
       case 2
       x2p = x;
       x2p(i) = x2p(i) + 2*h;
       g(i)= (-func(x2p) + 4*func(xp) -3*f)/(2*h);
       case 3
       xm = x;
       xm(i) = xm(i) - h;
       g(i)= (func(xp) - func(xm))/(2*h);
       
   end
           
end

end
