clc, clear

%% 迭代方法计算Legendre多项式
syms x
P0 = 1;
P1 = x;
fprintf('P0=%d, P1=%s\n',P0,P1);

xx=-1:0.01:1;

hold on
plot(xx,subs(P0,xx),'Linewidth',2)
plot(xx,subs(P1,xx),'Linewidth',2)

for n=1:4
    
    P = (2*n+1)/(n+1)*x*P1-n/(n+1)*P0;
    P0 = P1;
    P1 = P;
    plot(xx,subs(P1,xx),'Linewidth',2);
    fprintf('P%d=%s\n',n+1,simplify(P));
    
end

legend('P0','P1','P2','P3','P4','P5','Location','southeast')


%% 非迭代方法计算Legendre多项式

syms x
m = 5; % 生成次数小于m的所有Legendre多项式
k = 0:m;

M = k'*ones(1,m+1);

K = tril(ones(m+1,1)*k,0);

P_coeff = factorial(M+K)./(factorial(K)).^2./factorial(M-K); % 第k行为k次Legendre多项式的每项系数

P_vec = tril(P_coeff,0).*(x/2-1/2).^K; % 第k行为k次Legendre多项式的每项

P = simplify( sum(P_vec,2) );

disp(P);