eps = logspace(-1,-19);
x = 1.5;

% f(x) = x^(9/2)
n = 9/2;
df_ln = ((x+eps).^n - x.^n) ./ eps;
df_cp = imag((x+1i*eps) .^ n) ./ eps;

% real df(x) = 9/2 * x^(7/2)
realdf = n * x^(n-1);

loglog(eps, abs(df_ln - realdf), 'ko--');
hold on;
loglog(eps, abs(df_cp - realdf), 'bo-');

legend('finite difference', 'complex variable');
xlabel('epsilon');
ylabel('error in df(x)');
title('function f(x) = x^{9/2}');

l = loglog([1e-7, 1e-1], [1e-7, 1e-1].^2 * 2, 'r');
set(l, 'linewidth', 3);
text(1e-3/2, 1e-8, 'k = 2');

l = loglog([1e-8, 1e-1], [1e-8, 1e-1]*10, 'r');
set(l, 'linewidth', 3);
text(1e-6, 1e-6, 'k = 1');

l = loglog([1e-8/2, 1e-16], 0.3e-16./[1e-8/2, 1e-16]*10, 'r');
set(l, 'linewidth', 3);
text(0.3e-13, 1e-4, 'k = -1');
