%% Example 2.1
% Accuracy of Finite Difference Schemes

%% Define functions

% Function
f = @(x) sin(x)/x^3;
% Analytic first derivative
df = @(x) cos(x)/x^3-3*sin(x)/x^4;
% Forward difference
forward = @(x,h) (f(x+h) - f(x))/h;
% Central difference
central = @(x,h) ( f(x+h) - f(x-h))/2/h;
% Fourth order difference
fourth = @(x,h) ( f(x-2*h) - 8*f(x-h) + 8*f(x+h) - f(x+2*h) )/12/h;

%% Evaluate error in finite differences
x = 4;

H_for = logspace(-5,0,20);
H_cen = logspace(-5,0,10);
H_fth = logspace(-2.5,0,10);
for_err = [];
cen_err = [];
fth_err = [];

for i = 1:length(H_for)
    for_err(i) = abs(forward(x,H_for(i)) - df(x));
end

for i = 1:length(H_cen)
    cen_err(i) = abs(central(x,H_cen(i)) - df(x));
end

for i = 1:length(H_fth)
    fth_err(i) = abs(fourth(x,H_fth(i)) - df(x));
end

%% Generate plot
figure(1),clf
loglog(H_for,for_err,'o-')
hold on
loglog(H_cen,cen_err,'ko-')
loglog(H_fth,fth_err,'ro-')
xlabel('h -- grid spacing','FontSize',12)
ylabel('error','FontSize',12)
legend('1st order','2nd order','3rd order','Location','NorthWest')

%% Compute slopes
first_order = mean(diff(log(for_err))./diff(log(H_for)))
second_order = mean(diff(log(cen_err))./diff(log(H_cen)))
fourth_order = mean(diff(log(fth_err))./diff(log(H_fth)))
