clear

% cubic example for 4 points
x=(1:4)';
y=[3 -2 -5 0]';
A=ones(4,4);
A(:,2)=x;
A(:,3)=x.^2;
A(:,4)=x.^3;
a=A\y

t=(0:.1:5)';
p=polyval(flipud(a),t);
figure;plot(t,p,'k',x,y,'ro')

% Runge function example
N=15;
x=linspace(-1,1,N)';
y=1./(1+25*x.^2);
A=ones(N,N);
for k=2:N
    A(:,k)=x.^(k-1);
end
a=A\y;
t=(-1:.01:1)';
p=polyval(flipud(a),t);
figure;plot(t,p,'k',x,y,'ro',t,1./(1+25*t.^2),'b')
figure;plot(t,1./(1+25*t.^2)-p,'k',[-1 1],[0 0],'b')

% exp(x/2) example on [0,10] with nodes at x=0,2,...,10
x=(0:2:10)';
y=exp(x/2);
A=ones(6,6);
for k=2:6
    A(:,k)=x.^(k-1);
end
a=A\y;
t=(0:.1:10)';
pI=polyval(flipud(a),t);
figure;plot(t,pI,'b',t,exp(t/2),'k')
set(gca,'FontSize',20)
legend('Interpolating poly','Function exp(x/2)')

figure;plot(t,abs(exp(t/2)-pI),'b',t,abs(exp(5)/2^6/factorial(6)*(t-0).*(t-2).*(t-4).*(t-6).*(t-8).*(t-10)),'k')
set(gca,'FontSize',20)
legend('True error for interp poly','Theoretical upper bound')

pT=exp(5/2)*(1+(t-5)/2+(t-5).^2/2^2/factorial(2)+(t-5).^3/2^3/factorial(3)+...
    (t-5).^4/2^4/factorial(4)+(t-5).^5/2^5/factorial(5));
figure;plot(t,pT,'b',t,exp(t/2),'k')
set(gca,'FontSize',20)
legend('Taylor poly','Function exp(x/2)')

figure;plot(t,abs(exp(t/2)-pT),'b',t,abs(exp(5)/2^6/factorial(6)*(t-5).^6),'k')
set(gca,'FontSize',20)
legend('True error for Taylor poly','Theoretical upper bound')