%%
%% Listing 5.7 w książce "Obliczenia inżynierskie i naukowe", PWN 2011
%% (C) Piotr Krzyżanowski, 2011
%%
%% Możesz używać tego kodu zgodnie z licencją Creative Commons BY-SA, 
%% zob. http://creativecommons.org/licenses/by-sa/3.0/pl/
%%

function df = divdif(f, x, h)
% 
% DF = divdif(F, X, H)
%
% Wyznacza przybliżenie F'(X) ze wzoru (F(X+H)-F(X-H)) / (2H)
% Jeśli nie podano H, przyjmuje się H = max(|X|,1)*sqrt(eps).

if nargin < 3
	h = max(abs(x),1)*sqrt(eps);
end
f1 = feval(f,x+h);
f2 = feval(f,x-h);
df = (f1-f2)./(2*h);
end
