Function interpolation using Newton polynomial with Finite Differences

NewtonInterpolationFD

clear;
close all;
clc;

a = -10;
b = 10;
n = 14;

j = a:0.2:b;

h = (b - a)/n;

l = a:h:b;

syms x;
syms c;

f = sin(c);

y = subs(f, l);

P = y(1);

fdc(1:size(l, 2), 1:size(l, 2)) = 0;

for I=1:1:size(l, 2)
    for J=1:1:size(l, 2)-1
        if I == 1
            fdc(J, I) = y(J+1) - y(J);
        else
            fdc(J, I) = fdc(J+1, I-1) - fdc(J, I-1);
        end;
    end;
end;

qs = (x - l(1)) / h;
Q = qs;

for I=2:1:size(l, 2)
    P = P + fdc(1, I-1) * Q;
    Q = Q * (qs - I + 1) / I;
end;

plot(j, subs(f, j), '-b');
hold on;

plot(j, subs(P, j), '-r');
title('Num exp, n=7');
legend('f(x)', 'Interpolation');