n

%% Central Difference
% The Higher derivative Up to n " *Central Difference* " metrhod is implemented as
% shown here using MATLAB(R).

%% Clear Memory and Screen
clear
clc

%% Here, The user enter the order of derivative that he needs.
disp('Please Enter The order of Derivative');
Order = input('  Order = ');
disp('-------------------------------------');

%% Enter The Function.
fs = input('Please Enter The Function \n f(x)= ','s');
f  = inline(fs);

%% Input Mode
% Here, the user has to define where he wants to get the values of the 1st derivative.
% The user can input that using two formats: points that he wants to get the 1st derivative
% or the boundaries(a,b) of an period[a,b] besides a number of periods(n).
disp('-------------------------------------');
while 1
    disp('Please Choose points Input mode (Points or Boundary)');
    disp('Points like [X0 X1 X2 ... ]');
    disp('Boundary means enter [a,b] and n');
    Point_Mode = input('','s');
    disp('-------------------------------------');
    if ~(strcmpi( Point_Mode , 'Boundary' )||strcmpi( Point_Mode , 'Points' ))
        disp('Error!');
    else
        break;
    end;
end;
% Here, the code deals with the one of the two formats (Go back to Line 19).
if strcmpi( Point_Mode , 'Boundary' )
    a = input('Please Enter [a,b]\n  a = ');
    b = input('  b = ');
    n = input('Please Enter number of periods\n  n = ');
    h = (abs(b-a))/n;
    X = a:h:b;
    h(1:n+1) = h(1);
else
    X = input('Please Enter X \n For Example X = [0,1,5,6]\n  X = ');
    %Enter The Value of Step Size h.
    disp('Enter a value of (h) OR the corresponding (h) for every (x)');
    disp('Note: if numbers of (h) is less than number of (x), the remaining');
    disp('corresponding values of (h) will be last (h) entered');
    disp('-------------------------------------');
    h = input('Please Enter the value of Step Size\n  h = ');
    n = numel(X);
    nh = numel(h);
    % If data not complate for unequal we let the other h equal last value of h.
    if nh<n
        h(nh:n) = h(nh);
    end;
end;
disp('-------------------------------------');

% Find the Needed Derivative.
n = numel(X);
HD = X;
for j = 1:n
    HD(j) = 0;
    for i = 0:Order
        HD(j) = HD(j) + ((-1)^i) * nchoosek(Order,i) * feval(f,X(j)+(((Order/2)-i)*h(j)));
    end;
    HD(j) = HD(j)/(h(j)^Order);
end;

%Display X, h,  HD
fprintf('f(x) = %s\n',fs);
disp   ('+---------------+------------+-------------------+');
fprintf('|      x        |      h     |    %1.0f derivative   |\n',Order);
disp   ('+---------------+------------+-------------------+');
for i = 1:n
    fprintf('| %12.4f  |  %8.4f  |   %12.4f    |\n',X(i),h(i),HD(i));
end;
disp   ('+---------------+------------+-------------------+');