f = @(x)sin(pi*x).*exp(x);

% interpolation 
n = 4;
xx = -1:0.01:2;
% compute iterpolant
nodes =[-1,0,1,1.5,2];
coeffs = polyfit(nodes,f(nodes),n);
p3 = polyval(coeffs,xx);

plot(xx,f(xx),'k',nodes,f(nodes),'ro',xx, p3,'b', 'Linewidth', 2);
h_leg = legend('f(x)', 'interp. p3(x)');
set(h_leg,'FontSize',16); axis([-1,2,-5,2]);
pause

% compute optimal l2 interpolation with different weights
n = 4;
w0 = @(x)1+0.*x;
w1 = @(x)1+x.^4;
w2 = @(x)1./(1.2+x);
A = zeros(n+1,n+1);
b = zeros(n+1,1);
for i=0:n
    for j=0:n
    A(i+1,j+1) = integral(@(x)w0(x).*x.^(i+j),-1,2);
    end
    b(i+1) = integral(@(x)w0(x).*f(x).*x.^i,-1,2);
end
coeffs0 = A\b;
q0 = polyval(fliplr(coeffs0'),xx);
figure(1);
plot(xx,f(xx),'k', xx, q0,'m', 'Linewidth',2)
h_leg = legend('f(x)','optimal l2 with w=1');
set(h_leg,'FontSize',16);axis([-1,2,-5,2])
figure(2);
plot(xx,w0(xx),'r','Linewidt',2);
title('weight function')
pause

for i=0:n
    for j=0:n
    A(i+1,j+1) = integral(@(x)w1(x).*x.^(i+j),-1,2);
    end
    b(i+1) = integral(@(x)w1(x).*f(x).*x.^i,-1,2);
end
coeffs1 = A\b;
q1 = polyval(fliplr(coeffs1'),xx);
figure(1);
plot(xx,f(xx),'k', xx, q1,'m-.', 'Linewidth',2)
h_leg = legend('f(x)','optimal l2 with w focused around x=2');
set(h_leg,'FontSize',16);axis([-1,2,-5,2])
figure(2);
plot(xx,w1(xx),'r','Linewidt',2);
title('weight function')
pause

for i=0:n
    for j=0:n
    A(i+1,j+1) = integral(@(x)w2(x).*x.^(i+j),-1,2);
    end
    b(i+1) = integral(@(x)w2(x).*f(x).*x.^i,-1,2);
end
coeffs1 = A\b;
q2 = polyval(fliplr(coeffs1'),xx);
figure(1);
plot(xx,f(xx),'k',xx, q2,'m:', 'Linewidth',2);
h_leg = legend('f(x)','optimal l2 with w focused around x=-1');
set(h_leg,'FontSize',16);axis([-1,2,-5,2])
figure(2);
plot(xx,w2(xx),'r','Linewidt',2);
title('weight function')
pause


