Legrange's Interpolation

clc;
clear all;
syms p;
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);
% Divided difference formula
for i = 1:n
D(i,1) = y(i);
end
for i = 2:n
for j = 2:i
D(i,j)=(D(i,j-1)-D(i-1,j-1))/(x(i)-x(i-j+1));
end
end
fx = @(p) D(n,n);
for i = n-1:-1:1
fx = fx*(p-x(i)) + D(i,i);
end
fx = simplifyFraction(fx);
fprintf('Newton's iterated value using divided difference: f(%f) = %f \n', x0, subs(fx, p, x0));