clear;
format long;
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 diference 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);


%first order derivative
syms p;
product_p = 1;
sum = 0;
for i=1:range-x
    product_p = product_p * (p-(i-1));
    diff_p = diff(product_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(first);