%**************************************************************************
%*     Problem 2.2                                                        *
%*                                                                        *
%*  This program evaluates the second-order derivative of a given         *
%*  function f(x) at x=1 using the 3-point and 5-point methods.           *
%*  The derivative and numerical errors are plotted.                      *
%*                                                                        *
%*     Programed by Ryoichi Kawai for Computational Physics Course.       *
%*     Last modification:  01/31/2016.                                    *
%**************************************************************************
clear all;

% define the function
f = @(x) x^4/12;

x=1;  % the point at which derivative is evaluated.
imax=10; % number of different displacements

% title and headers
display('Absolute Errors')
fprintf('%4s %18s %18s\n','h','3-points','5-points')

for i=1:imax
    
    % Small displacement
    h(i)=10^(-i+1);
    
    % Evaluation of numerical derivative
    s3(i)=(f(x+h(i))+f(x-h(i))-2*f(x))/h(i)^2;
    s5(i)=(16*(f(x+h(i))+f(x-h(i)))-f(x+2*h(i))-f(x-2*h(i))-30*f(x))/(12*h(i)^2);
   
    % Errors
    err_s3(i)=abs(1-s3(i));
    err_s5(i)=abs(1-s5(i));
    
    % Display the errors
    fprintf('%6.1e %18.10e %18.10e\n',...
            h(i),s3(i),s5(i));
end

h2 = h.^2; % h^2 (the error of the mean value formula is order of h^2)
h4 = h.^4;

% Plot data
subplot(1,2,1)  % left panel
q=semilogx(h(1:imax),s3(1:imax),'-o',...
    h(1:imax),s3(1:imax),'-s');
title('Numerical Derivative');
xlabel('h','fontsize',14);
ylabel('Derivative','fontsize',14);
set(q(1),'Color','blue');
set(q(2),'Color','red');
legend(q,{'3-point','5-point'});
legend(q,'Location','NorthWest');

subplot(1,2,2) % right panel
p=loglog(h(1:imax),err_s3(1:imax),'o',...
    h(1:imax),err_s5(1:imax),'s',...
    h(1:imax/2),h2(1:imax/2),'--',h(1:imax/2),h4(1:imax/2),'--');
title('Errors in the finite difference methods');
xlabel('h','fontsize',14);
ylabel('absolute error','fontsize',14);
set(p(1),'Color','blue');
set(p(2),'Color','red');
set(p(3),'Color','blue','LineWidth',1);
set(p(4),'Color','red','LineWidth',1);
legend(p,{'3p','5p','h^2','h^4'});
legend(p,'Location','SouthWest');
hold off