%interpolates sin(x) at n equally spaced points in [0,1]. Also evaluates 
%sin(x) and interpolating polynomial P(x) at 1000 equally spaced points in [0,1]

f = @(x) sin(x);
    
n = 3;
x = linspace(0,1,n);    %equally spaced nodes in [0,1]
c = zeros(1,n);         %coefficients of Newton interpolant

for i=1:n
    c(i) = f(x(i));
end

for j=2:n
    for i=n:-1:j
        c(i) = (c(i)-c(i-1))/(x(i)-x(i-j+1));
    end
end

x_val = linspace(0,1,1000);  %interpolant P is evaluated on x_val
y_int = zeros(1,1000);

for i=1:1000
    y_int(i) = c(1);    %inner loop computes y_int(i) = P(x_val(i))
    p = 1;
    for j=2:n
        p = p*(x_val(i) - x(j-1));
        y_int(i) = y_int(i) + c(j)*p;
    end
end                           %y_int = interpolated values of f on x_val

y_val = zeros(1,1000); 

for i=1:1000
    y_val(i) = f(x_val(i));   %y_val = exact values of f on x_val
end                           

disp('the actual maximum error is...')
max(abs(y_int-y_val))

disp('the theoretical error bound is...')
(4*n)^(-1)*(1/(n-1))^n
