function dlr = demo_lagrange_runge(degree,excent)
%
% Interpolation of the Runge function by Lagrange polynomial and splines
%    
% degree = No. of interp. intervals (= polynomial degree)
% excent = Center of interpolation shifted w.r.t. symmetry center of function
%    
    
% Function to be approximated    
f=inline('[1./(1+25*x.*x)]','x');

% No. of points for plotting
nd = 200;

% No. of interpolation intervals
nomg =  degree; 
if (excent==0)
    % Range for interpolation and plotting
    xap_min = -1.0;  xap_max = 1.0;
    xpl_min = -1.1;  xpl_max = 1.1;
else
    xap_min = -1.5;  xap_max = 0.5;
    xpl_min = -1.6;  xpl_max = 0.6;
end    
ypl_min = -0.4;  ypl_max = 1.1;
loc = 'NW';
  
% Points for plotting
dxpl = (xap_max-xap_min) / nd;
xpl  = xpl_min:dxpl:xpl_max;

% Interpolation points for Lagrange and spline
domg = (xap_max-xap_min) / nomg;
xi = xap_min:domg:xap_max;

% Values of polynomial omega at each plotting point
% version 1
omega = xpl - xi(1);
for k=2:nomg+1
    omega = omega.*(xpl-xi(k));
end
% version 2
%amat = diag(x);
%pol0 = poly(amat);
%omega = polyval(pol0,xpl);


figure

subplot(2,1,1)

% plot values of omega
plot(xpl,omega);
% plot axes
yop_min=min(omega(nd/20:19*nd/20));
yop_max=max(omega(nd/20:19*nd/20));
hold on;
plot([xpl_min, xpl_max],[0 0],'k--');
plot([0 0],[yop_min yop_max],'k--');
hold off;
axis([xpl_min xpl_max yop_min yop_max]);
title('Polynomial omega(x)');

% Analyt. values: "continuous" (for plotting) and discrete (x_i,y_i)
yi  = f(xi);
ypl = f(xpl);
% Interpolation by Lagrange and spline
p  = polyfit(xi,yi,nomg);
pg = polyval(p,xpl);
ys = spline(xi,yi,xpl);

subplot(2,1,2)
% plot the interpolations
plot(xpl,ypl, 'k', xi,yi, 'k+', xpl,ys, 'r--', xpl,pg, 'b')
legend({'analytic';'discrete y_i';'spline';'Lagr. pol.'},'location',loc)
title('Interpolation of y = 1/(1+25 x^2)')
axis([xpl_min xpl_max ypl_min ypl_max])
% plot axes
hold on;
plot([xpl_min, xpl_max],[0 0],'k--');
plot([0 0],[ypl_min ypl_max],'k--');
hold off;

