clc,clear

%% 数据（data）
x0 = -5:5;
y0 = 1./(1+x0.^2);

x0_Chebyshev = 5*cos( (2*(0:10)+1)*pi/22 );
y0_Chebyshev = 1./(1+x0_Chebyshev.^2);

Ln = Lag(x0,y0);
Ln_Chebyshev = Lag(x0_Chebyshev,y0_Chebyshev);


%% 画图（plot figure）
hold on;
plot(x0,y0,'ro','Linewidth',2);
plot(-5:0.01:5,1./(1+(-5:0.01:5).^2),'r','Linewidth',2);
plot(-5:0.01:5,subs(Ln,-5:0.01:5),'b--','Linewidth',1);
plot(-5:0.01:5,subs(Ln_Chebyshev,-5:0.01:5),'k-','Linewidth',1);

legend('data','f(x)','Lagrange','Lagrange-Chebyshev','Location','best')


function Ln = Lag(x0,y0)
%% 构造Lagrange函数（Get Lagrange polynomial）
    syms x

    omega_vec = x*ones(size(x0)) - x0;

    omega = prod(omega_vec);

    omega_der = diff(omega);

    Ln_vec = y0.*(omega./(x-x0))./subs(omega_der,x0);

    Ln = sum(Ln_vec);
    
end