1-6

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

%% Clear Memory and Screen
clear
clc

%% Firstly, the user has to choose between two input modes
% 1- Function: Formula of a function.
% 2- Data    : some (x,y) points.
% the user has to write one of two words (Function or Data) to indicate the mode.
while 1
    disp('Please Choose your Input mode (Function or Data)');
    Mode = input('','s');
    disp('-------------------------------------');

    if ~(strcmpi( Mode , 'Function' )||strcmpi( Mode , 'Data' ))
        disp('Error!')
    else
        break;
    end;
end;
%% Enter the order of Derivative
while 1
    disp('Please Enter The order of Derivative (1 to 6)');
    Order = input('  Order = ');
    disp('-------------------------------------');
    if (Order ~= 1)&&(Order ~= 2)&&(Order ~= 3)&&(Order ~= 4)&&(Order ~= 5)&&(Order ~= 6)
        disp('Error!');
    else
        break;
    end;
end;

%% Function Mode
% The user has chosen the Function mode (Go to line 10 ).
if strcmpi( Mode , 'Function' )
    % Enter the function
    fs = input('Please Enter The Function \n f(x)= ','s');
    f  = inline(fs);
    disp('-------------------------------------');

    % 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).
    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 44).
    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 Needed Derivative.
    n = numel(X);
    HD = X;
    for i = 1:n
        if Order == 1
            HD(i) = (feval(f,X(i)+h(i))-feval(f,X(i)-h(i)))/(2*h(i));
        elseif Order == 2
            HD(i) = (feval(f,X(i)+h(i))+feval(f,X(i)-h(i))-2*feval(f,X(i)))/(h(i)*h(i));
        elseif Order == 3
            HD(i) = (feval(f,X(i)+2*h(i))-2*feval(f,X(i)+h(i))+2*feval(f,X(i)-h(i))-feval(f,X(i)-2*h(i)))/(2*h(i)^3);
        elseif Order == 4
            HD(i) = (feval(f,X(i)+2*h(i))-4*feval(f,X(i)+h(i))+6*feval(f,X(i))-4*feval(f,X(i)-h(i))+feval(f,X(i)-2*h(i)))/(h(i)^4);
        elseif Order == 5
            HD(i) = (feval(f,X(i)+3*h(i))-4*feval(f,X(i)+2*h(i))+5*feval(f,X(i)+h(i))-5*feval(f,X(i)-h(i))+4*feval(f,X(i)-2*h(i))-feval(f,X(i)-3*h(i)))/(2*h(i)^5);
        elseif Order == 6
            HD(i) = (feval(f,X(i)+3*h(i))-6*feval(f,X(i)+2*h(i))+15*feval(f,X(i)+h(i))-20*feval(f,X(i))+15*feval(f,X(i)-h(i))-6*feval(f,X(i)-2*h(i))+feval(f,X(i)-3*h(i)))/(h(i)^6);
        end;
    end;

    %Display X, h, Needed Derivative.
    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   ('+---------------+------------+-------------------+');
end;

%% Data Mode
% The user has chose the Data mode (Go back to line 11 ).
if strcmpi( Mode , 'Data' )
    % The user has to enter the Boundary of interval.
    a = input('Please Enter the interval [a,b]\n  a = ');
    b = input('  b = ');
    n = input('Please Enter number of periods\n  n = ');
    if b<a
        tmp = a;
        a = b;
        b = tmp;
    end;
    h = (b-a)/n;
    X = a:h:b;

    % The user has to enter Y(x).
    Y = X;
    for i = 1:n+1
       fprintf('Please Enter Y(%1.4f) = ',X(i));
       Y(i) = input('');
    end;

    % Find Needed Derivative.
    n = numel(X);
    HD = X;
    HP = floor((Order+1)/2); % half of points that we need to find the drivative
    for i = HP+1 : n-HP
        if Order == 1
            HD(i) = (Y(i+1)-Y(i-1))/(2*h);
        elseif Order == 2
            HD(i) = (Y(i+1)+Y(i-1)-2*Y(i))/(h*h);
        elseif Order == 3
            HD(i) = (Y(i+2)-2*Y(i+1)+2*Y(i-1)-Y(i-2))/(2*h^3);
        elseif Order == 4
            HD(i) = (Y(i+2)-4*Y(i+1)+6*Y(i)-4*Y(i-1)+Y(i-2))/(h^4);
        elseif Order == 5
            HD(i) = (Y(i+3)-4*Y(i+2)+5*Y(i+1)-5*Y(i-1)+4*Y(i-2)-Y(i-3))/(2*h^5);
        elseif Order == 6
            HD(i) = (Y(i+3)-6*Y(i+2)+15*Y(i+1)-20*Y(i)+15*Y(i-1)-6*Y(i-2)+Y(i-3))/(h^6);
        end;
    end;

    %Display X, Y, H, HD
    disp   ('+---------------+---------------+---------------------+');
    fprintf('|      x        |      y        |     %1.0f derivative    |\n',Order);
    disp   ('+---------------+---------------+---------------------+');
    for i = 1 : n
        if (i<HP+1)||(i>n-HP)
            fprintf('| %12.4f  | %12.4f  | Can`t be calculated |\n',X(i),Y(i));
        else
            fprintf('| %12.4f  | %12.4f  |     %12.4f    |\n',X(i),Y(i),HD(i));
        end;
    end;
    disp   ('+---------------+---------------+---------------------+');

end;