legendpoly.m

to caluculate the legendre polynomials

% Given nonnegative integer n, compute the 
% Legendre polynomial P_n. Return the result as a vector whose mth
% element is the coefficient of x^(n+1-m).


function lp = LegendPoly(n)

if n==0 
    lp = 1;
elseif n==1
    lp = [1 0]';
else
    
    lpm2 = zeros(n+1,1);
    lpm2(n+1) = 1;
    lpm1 = zeros(n+1,1);
    lpm1(n) = 1;

    for k=2:n
        
        lp = zeros(n+1,1);

        for e=n-k+1:2:n
            lp(e) = (2*k-1)*lpm1(e+1) - (k-1)*lpm2(e);
        end
        
        lp(n+1) = lp(n+1) - (k-1)*lpm2(n+1);
        lp = lp/k;
        
        if k<n
            lpm2 = lpm1;
            lpm1 = lp;
        end
        
    end
    
end