%**************************************************************************
%*     Example  2.2                                                       *
%*     filename: ch02pr02.m                                               *
%*     program listing number: 2.2                                        *
%*                                                                        *
%*  This program evaluates the derivative of a given function func(x)     *
%*  at x=1 using the mean finite difference method with the accuracy      *
%*  specified by tolerance.                                               *
%*                                                                        *
%*     Programed by Ryoichi Kawai for Computational Physics Course.       *
%*     Last modification:  01/11/2015.                                    *
%**************************************************************************
clear all;

% Set the tolerance
tol=0.0001  % 

% Set the number of points to evaluate the derivative
N=101;
% Find the grid size
dx=2*pi/(N-1);
% Set the points where the derivative is evaluated.
x = (0:N-1)*dx;

% initial 
h=2; % initial interval

% any value bigger than tol is OK here.
err = tol+1;

fprintf('%5s %18s\n','h','max error')

% Repeat until error is smaller than tolerance.
while err>tol
    % set the interval
    h=h/2;
    % Derivative by the mean finite difference method
    numerical=(sin(x+h)-sin(x-h))/(2*h);
    % Exact value
	exact=cos(x);
    % Largest absolute error
    err=max(abs(numerical-exact));
    % Show the results.
    fprintf('%8.6f  %15.6e\n',h,err)
end
fprintf('Tolerance is OK. Good value of h=%10.6f\n',h)

% Plot the results.
subplot(1,2,1)
p=plot(x,numerical,'o',x,exact);
set(p(1),'Linewidth',2,'Color','red')
set(p(2),'Linewidth',2,'Color','blue')
xlabel('x','Fontsize',14)
ylabel('$\frac{d}{dx} \sin(x)$','Interpreter','LaTeX','Fontsize',14) 
legend('numerical','exact')
legend('Location','Southeast')

subplot(1,2,2)
q=plot(x,abs(numerical-exact));
set(q(1),'Linewidth',2,'Color','black')
xlabel('x','Fontsize',14)
ylabel('Absolute error','Fontsize',14) 

