Produttoria di funzioni su matlab

[maledetto]

salve a tutti, sono alle prese con una relazione di calcolo numerico. L'argomento è la determinazione dei coefficienti di newton-cotes per la risoluzione degli integrali. Siccome per determinare i coefficienti devo risolvere un integrale implemento nel mio script anche il metodo di simpson. Nella prima parte dello script dovrei fare la produttoria di funzioni ma non so proprio come fare. Ho provato a prendere un vettore ones ma credo che non si possa sostituire un coefficiente di un vettore con una funzione. Ecco ciò che ho fatto
function[ckn] = newtcotes(b,m)
a=0;
h=(b-a)/(2*m);
for k=0:b
sum=0;
ckn1=ones(1,b+2);
for j=0:b
if (j~=k)
ckn1(1,j+2)=@(x)(((x-j)/(k-j)).*ckn1(1,j+1));
end
end