function d2p=poly_d2l(m,c)
%% values of second derivative of lagrange polynomials through 0:1/m:1 at c
% function d2p=poly_d2l(m,c);
% INPUT:
%       m lagrange degree
%       c evaluation point(s) in [0,1]
% OUTPUT:
%       d2p values of second derivative of lagrange polynomials 
%           through 0:1/m:1 at c

% (c) DDE-BIFTOOL v. 3.1.1(19), 11/04/2014 
%
% 
%
d2p=zeros(m+1,length(c));
if m==2
    d2p(1,:)=4;
    d2p(2,:)=-8;
    d2p(3,:)=4;
elseif m==3
    d2p(1,:)=-13.5*(c-0.8591167563965418+c-0.4742165769367915);
    d2p(2,:)=40.5*(c-0.8495279234516212+c-0.2615831876594899);
    d2p(3,:)=-40.5*(c-0.7384168123405100+c-0.1504720765483788);
    d2p(4,:)=13.5*(c-0.5257834230632086+c-0.1408832436034581);
elseif m==4
    d2p(1,:)=42.666666666666667*((c-0.625).*(c-0.9045084971874737)+ ...
        (c-0.3454915028125263).*(c-0.9045084971874737)+ ...
        (c-0.3454915028125263).*(c-0.625));
    d2p(2,:)=-170.666666666666667*((c-0.61727336112119537025).* ...
        (c-0.90181300471826281827)+ ...
        (c-0.16841363416054181147).*(c-0.90181300471826281827)+ ...
        (c-0.16841363416054181147).*(c-0.61727336112119537025));
    d2p(3,:)=256*((c-0.5).*(c-0.8952847075210475)+ ...
        (c-0.1047152924789526).*(c-0.8952847075210475)+ ...
        (c-0.1047152924789526).*(c-0.5));
    d2p(4,:)=-170.666666666666667*((c-0.38272663887880462974).* ...
        (c-0.83158636583945818852)+ ...
        (c-0.09818699528173718172).*(c-0.83158636583945818852)+ ...
        (c-0.09818699528173718172).*(c-0.38272663887880462974));
    d2p(5,:)=42.666666666666667*((c-0.375).*(c-0.6545084971874737)+ ...
        (c-0.09549150281252627).*(c-0.6545084971874737)+ ...
        (c-0.09549150281252627).*(c-0.375));
else
    t=0:1/m:1;
    for j=1:m+1
        for k=1:m+1
            if k~=j
                ss=zeros(1,length(c));
                for s=1:m+1
                    if s~=j && s~=k
                        f=ones(1,length(c));
                        for l=1:m+1
                            if l~=k && l~=j && l~=s
                                f=f.*(c-t(l))/(t(j)-t(l));
                            end
                        end
                        ss=ss+f/(t(j)-t(s));
                    end
                end
                d2p(j,:)=d2p(j,:)+ss/(t(j)-t(k));
            end
        end
    end
end
end
