% Generates Lagrangian Interpolant and Spline Fit to m Fit Points, Error for both, and Comparison of Error
% David Bloore, July 2011

% first, let's get the user's input

m=input('Number of Mesh Points?');

% this generates the function by establishing a domain, then calculating f(x)

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

% again we generate the the fit points, first the values for x, then f(x)

xin=linspace(-1+2*(m+1)^-1,1-2*(m+1)^-1,m);
yin=(1+25*xin.^2).^-1;

% interpolant

pp=polyfit(xin,yin,m-1);
y1=polyval(pp,x);

% spline

yB=spline(xin,[0 yin 0]);
y2=ppval(yB,x);

% errors (computed as a simple difference)

interror=y1-f;
splinerror=y2-f;
errordifference=interror-splinerror;

% plot and parameters

plot(x,f,'-',xin,yin,'o',x,y1,'--',x,y2,'-.',x,interror,'--',x,splinerror,'-.',x,errordifference,':');
legend('f(x)=(1+25x^2)^-1','Fit Points','Interpolant of degree m-1','B-spline','Error of Interpolant','Error of B-spline','Error Int - Error B');
xlabel('-1 \leq x \leq 1');
ylabel('f(x)=(1+25x^2)^-1');
title('Interpolation and Spline Through f(x)=1/(1+25x^2) -- David Bloore, 18.085, July 28,2011');
set(gca,'XTick',-1:(m+1)^-1:1);
ylim([-0.25 1.1]);
grid on;