NumericalDifferentiationAugust2019

clear all;
clc;
format long;

%input data points
x=[1 2 3 4 5];
y=[2 5 9 12 20];
x_at = 1;

range=length(x);
h=x(2)-x(1);
p_val=(x_at-x(1))/h;

%finding index of x_at
x_index=1;
for i=1:range
    if(x_at==x(i))
        x_index = i;
        break;
    end
end


%forward difference table
new_y=y(x_index:range);
diff_y=zeros(1, range - x_index);
for i=1:range - x_index
    temp=diff(new_y, i);
    diff_y(i)=temp(1);
end

disp(diff_y);
syms p;
prod_p = 1;
sum=0;
%First Order derivative
for i=1:range - x_index
    prod_p = prod_p * (p-(i-1));
    diff_p = diff(prod_p, 1);
    co_eff = (diff_p/factorial(i))*diff_y(i);
    eachTermVal=vpa(subs(co_eff, p, p_val));
    sum=sum+eachTermVal;
end

first=sum/h;

disp(eachTermVal);