%This script will compute polynimial interpolants to e^x at 2^n equally
%spaced nodes for n= 1-5 and compare the polynomial to the values of e^x at
%501 equally spaced points, finding the maximum error at those points.  It
%produces a plot of each polynomial and a plot of the error vs. n.
x(1)=-1;
y(1)=exp(-1);
%These values are the same for all the polynomials so....It would of course
%be more efficient to produce the nodes in the order -1, 1, 1/2 14, 24,
%1/8, 3/8, 5/8, 7/8, ... then you could just call the DivDiff routine once
%and then read off the coref's to the different polynomials.  But the
%computer had much more time than I did so, I did not program it that
%efficiently.
for i=1:5
    n(i)=2^i;
    for j=1:n(i)
        x(j+1)=-1+2*j/n(i); %compute the nodes
        y(j+1)=exp(x(j+1)); %compute the y-values at the nodes
    end
    A=DivDiff(x,y,n(i)); %compute the divided dffierence coef's
    max(i)=0; %initialize teh max error to 0
    for j=0:499
        xx(j+1)=-1+j*.004; %get the next of 501 equally spaced points 
        %(actually you don't need the first and last since they are nodes).
        pval(j+1)=Interp(A(1,1:n(i)+1),x,n(i),xx(j+1)); %evaluate the polynomial
        error=abs(exp(xx(j+1))-pval(j+1)); %find the absolute error at this point
        if (error > max(i))
            max(i)=error; %update the max error if this point's is larger.
        end
    end
    figure %start a new plot window
    plot(xx,pval) %plot the interpolating polynomial
end
figure %start a new plot window
plot(n,max) %plot the error, None of the plots were that interesting....
format long
disp(n) %display the max error and n
disp(max)


    