% Generates a Series of Lines of Error of Interpolant Fits to m Points of 1/(1+25x^2)
% David Bloore, July 2011

% user-defined parameters

m=input('Lower bound of m?');
r=input('Length of interval of m?');
d=input('Number of discrete points along x axis?');

% blank matrices for plot3

M=linspace(m,m+r,r+1);
X=zeros(d,r+1);
Y=zeros(d,r+1);
Z=zeros(d,r+1);

% domain and range of interest

x=linspace(-1,1,d);
f=(1+25*x.^2).^-1;

% this loop computes the error of the interpolant for a specific value of m,
% assigns the appropriate x,y,z values to the appropriate columns of X,Y,Z

counter=1;

while counter<r+2;
    
    mtemp=m+counter-1;
    Mtemp=mtemp*ones(1,d);
    
    xin=linspace(-1+2*(mtemp+1)^-1,1-2*(mtemp+1)^-1,mtemp);
    zin=(1+25*xin.^2).^-1;
    
    pp=polyfit(xin,zin,mtemp-1);
    zBe=polyval(pp,x)-f;
    
    X(1:d,counter)=x(1:d);
    Y(1:d,counter)=Mtemp;
    Z(1:d,counter)=zBe(1:d);
    
    counter=counter+1;
end

% each column of Z corresponds to a line, same for X,Y--so the first column of each matrix produces the first line, and so on

plot3(X,Y,Z);
xlabel('-1 \leq x \leq 1');
ylabel('Number of Points Used for Interpolant Fit');
zlabel('Error = Interpolant fit - 1/(1+25x^2)');
title('Series of Lines of Error of Interpolant Fits to m Points of 1/(1+25x^2)');
set(gca,'XTick',-1:1/(m+1):1);
set(gca,'YTick',m:m+r);
set(gca,'LineStyleOrder', '--|-.|:');
set(gca,'YDir','reverse');
view(-60,25);
% zlim([-2 2]);
grid on;
hidden off;