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';