% approximate derivative
% with difference method
%

 
function y = Dapprox(func,u)



abu = eye(length(u));

h = sqrt(eps);

for i=1:1:length(u)
    
  y(:,i) = (func(u+h*abu(:,i)) -func(u - h*abu(:,i)) )/(2*h);
  
end
