matlab Newton-Cotes数值积分公式

function c=cotescoeff(n) %cotes系数
syms t
for i=1:n+1
    k=i-1;
    for j=1:n
       Ck(i)=(-1)^(n-k)/factorial(k)/factorial(n-k)/n*quad1(@(t)intfun(t,n,k),0,n)
    end
end

function [y,ck,ak]=newtoncotes(fun,a,b,n)   %牛顿-科特斯数值积分公式
% n=1,即梯形公式
%n=2，即辛普森公式
%n=4，即科特斯公式
xk=linspace(a,b,n+1);
fx=fun(ck);
ck=cotescoeff(n);  %科特斯系数
ak=(b-a)*ck;
y=ak*fx';