function yInt = NewtonInterpolation(x, y, xx)
% NewtonInterpolation(x, y, xx): Uses an (n-2)order Newton interpolating
% polynomial based on n data points (x,y)
% to determine a value of a dependent variable (yInt)
% at a given value of the independent variable, xx
% input:
% x - independent variable (row matrix)
% y - dependent variable (row matrix)
% xx - value of independent variable at which interpolation is calculated
% output:
% yInt - interpolated value of dependent variable

%compute
n = length(x);
if length(y) ~= n
    error('x and y must have the same length');
end
b = zeros(n,n);
%assign dependent variables to the first column of b
b(:,1)=y(:); %ensures that y is a column vector
for j = 2:n
    for i = 1:n-j+1
        b(i,j) = (b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i));
    end
end

fprintf('\n\t\t\t\t\t\t\tx\t\t\t\t\t\t\tfx');
for k = 1:n-1

    if k==1
        fprintf('%30.fst',k);
    end
    if k==2
       fprintf('%30.fnd',k);
    end
    if k==3
        fprintf('%30.frd',k);
    end
    if k>3
        fprintf('%30.fth',k);
    end

end
l = b';
v = [x;l];
result = v';
for k = 1:n
    fprintf('\n');
    fprintf('%30.5f\t',result(k,:));
end

%uses the finite divided difference to interpolate
xt = 1;
yInt = b(1,1);

fprintf('\n\n\ty =');
fprintf('% 10.5f\t',b(1,1));

for j = 1:n-1
    xt = xt*(xx-x(j));
    yInt = yInt + b(1,j+1)*xt;
        fprintf('+');
        fprintf('% 10.5f*(x-%f)',b(1,j+1),x(j));
    if j == 1
        fprintf('\n\t\t\t\t\t');
    end
    if j>1
        k = 1;
        p=j;
        while k<p
            fprintf('*(x-%f)',x(p-1));
            p = p-1;
        end
        if j < n-1
            fprintf('\n\t\t\t\t\t');
        end
    end
end