function fd = complex_step(f,x,h)
%COMPLEX_STEP  Complex step approximation to derivative.
%  fd = COMPLEX_STEP(f,x,h) computes the complex step approximation 
%  fd to the derivative of f at x, using step h (default 1e-100).
if nargin < 3, h = 1e-100; end
fd = imag( f(x + sqrt(-1)*h) )/h;