Richardson Extrapolation

%% Richardson Extrapolation
% The first derivative " *Richardson Extrapolation* " method is implemented as
% shown here using MATLAB(R).

%% Clear Momory and Screen
clear
clc
%% Enter The Function
fs = input('Please Enter The Function \n f(x)= ','s');
f  = inline(fs);
disp('-------------------------------------');

%% Choose Output Mode
% 1- Direvative table (Dn1, Dn2, Dn3, ....) for only one point.
% 2- only the value of Derivative for many value of x.
while 1
    disp('choose output mode');
    disp('1 => Direvative table (Dn1, Dn2, Dn3, ....) for only one point');
    disp('2 => only the value of Derivative for many value of x');
    Out_Mode = input('Enter mode (1 or 2) = ');
    disp('-------------------------------------');
    if Out_Mode == 1 || Out_Mode == 2
        break;
    else
        disp('Error');
    end;
end;

%% 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).

% 1- Direvative Table Mode.
if Out_Mode == 1
    Point_Mode = 'Points'; %Just one Point X(1).
end;
% 2- Value Mode.
while Out_Mode == 2
    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;

%% Enter Boundary or Points :: Values of X
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 the value of x = ');
end;
disp('-------------------------------------');

%% Choose the accuracy (h,m)
disp('Choose the accuracy you need');
disp('If the T.E = O(h^m) and h<1 then for greater m we get small error');
while 1
    disp('Note: m must be even number');
    m = input('Enter m = ');
    disp('-------------------------------------');
    if (mod(m,2) == 0)
        break;
    end;
end;
m = int32(m/2);

%% Enter The Value of Step Size (h) for points mode
if strcmpi( Point_Mode , 'Points' )
    if Out_Mode == 2
        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('-------------------------------------');
    end;
    h = input('Please Enter the value of Step Size\n  h = ');
    disp('-------------------------------------');
    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;

%% Find First Drivative

D = zeros(m); % prepare size of the matrix D to avoid resizing.
% First Derivative of X(j) ==> FD
FD = X; % prepare size of the matrix FD to avoid resizing.

% Direvative Table Mode just use X(1).
% This Condition to deal with more than one value of X.
if Out_Mode == 1
    n = 1;
else
    n = numel(X);
end;

for j = 1:n
    hj = h(j);
    D(1,1) = (feval(f,X(j)+hj)-feval(f,X(j)-hj))/(2*hj);
    for i = 1:1:m-1
        hj = hj/2;
        D(i+1,1) = (feval(f,X(j)+hj)-feval(f,X(j)-hj))/(2*hj);
        for k = 1:1:i
            D(i+1,k+1) = D(i+1,k)+(D(i+1,k)-D(i,k))/((4.0^double(k))-1);
        end;
    end;
    FD(j) = D(m,m);
end;

%% Output for Table Mode
if Out_Mode == 1
    fprintf('f(x) = %s\n',fs);
    Line = '+------+------------+';
    for i = 1:m
        Line = strcat(Line,'-----------------+');
    end;

    disp(Line); % Before Header.
    fprintf('|   i  |      h     |');
    for i = 1:m
        fprintf('       Di,%1.0d      |',i);
    end;
    fprintf('\n');

    disp(Line); % Before Data.
    for i = 1:m
        fprintf('|  %2.0f  |  %8.4f  |',i,h);
        for j = 1:m
            if j>i
                fprintf('                 |');
            else
                fprintf('   %12.4f  |', D(i,j));
            end;
        end;
        fprintf('\n');
        h =h/2;
    end;

    disp(Line); % End.
end;

%% Output for Value Mode
if Out_Mode == 2
    fprintf('f(x) = %s\n',fs);
    disp   ('+---------------+------------+-------------------+');
    fprintf('|      x        |      h     | first derivative  |\n');
    disp   ('+---------------+------------+-------------------+');
    for i = 1:n
        fprintf('| %12.4f  |  %8.4f  |   %12.4f    |\n',X(i),h(i),FD(i));
    end;
    disp   ('+---------------+------------+-------------------+');
end;