Newton's Raphson Backward Difference

clc;
clear all;
syms p
%Defining x and y arrays
x = input('Input x values as an array : ');
y = input('Input y values as an array : ');
x0 = input('Enter the x for which y is to be found : ');
%Finding the length of x
n = length(x);
D = zeros(n,n);
D(:,1) = y;
for j = 2:n
for i = n:-1:j
D(i,j) = D(i,j-1)-D(i-1,j-1);
end
end
h = x(2) - x(1);
u = (p-x(n))/h;
fx = @(p) y(n);
G = u;
for k = 1:n-1
fx = fx+G*D(n,k+1);
G = G*(u+k)/(k+1);
end
fx = simplify(fx)
fprintf('Value of f(x0) using Newtons Backward interpolation formula is f(%f) = %f',x0,subs(fx,p,x0))