MatLab / Newton

function Newton(X, X1, f)

for i=1:1:length(X)
F(i)=subs(f,X(i));
end

A(1) = F(1);
A(2) = (F(2) - A(1))/(X(2)-X(1));
for i=3:1:length(F)
B = A(1);
for j = 2:1:i-1
C = 1;
for k = 1:1:j-1
C = C * (X(i) - X(k));
end
B = B + A(j)*C;
end
A(i) = (F(i) - B)/(C * (X(i) - X(i-1)));
end

syms x
p = A(1);
P(1) = A(1);
for i = 2:1:length(X)
C = 1;
for j = 1:1:i-1
C = C * (x - X(j));
end
p = p + A(i)*C;
end
P = subs(p, X1);
plot(X1, P)
P = vpa(P,4)

Y = 0:pi/100:1;
R = abs(subs(f, Y) - subs(p, Y));
A = max(R);
vpa(A)
end