function [ dy ] = Lyness( f, zeta, r, n )
%LYNESS Differentiate f by complex evaluation around zeta
%   Spectrally convergent formula from Lyness and Moler 1967
%   dy  = Lyness( f, zeta, r, n ) gives dy = f'(zeta) + O(r^n)
tau = zeta + r* exp(2*pi*1i*(0:n-1)/n);
y = feval( f, tau );
dy = sum(y.*exp(-2*pi*1i*(0:n-1)/n))/n/r;
end

