%% composite.m %% Compute the weights of the Newton-Cotes quadrature quadrature% on [0,1] with equidistant nodes.%%% errror functionf = inline('exp(-x.^2)*2/sqrt(pi)', 'x');a = 0; b = 2;int = erf(b);% bessel functionf = inline('cos( 8*sin(x) - x )/pi', 'x');a = 0; b = pi;int = bessel(1,8);% f = inline('sin(x).^2.*cos(x).^4', 'x');a = 0; b = 2*pi;int = b/16;% f = inline('4./(1+x.^2)', 'x');a = 0; b = 1;int = pi;% composite trapezoidal rulefor n = 1:16    % subinterval length    h = (b-a)/n;    % equidistant nodes    x = (a:h:b);    % function values at equidistant nodes    fx = feval( f, x);    % composite Trapezoidal rule    T(n) = h*( 0.5*(fx(1)+fx(n+1)) + sum(fx(2:n)) );end% composite Simpson rulefor n = 1:18    % subinterval length    h = (b-a)/n;    % equidistant nodes and midpoints    x  = (a:h:b);    % function values at equidistant nodes    fx = feval( f, x);    % composite Simpson rule    S(n) = (h/6)*( fx(1)+fx(n+1) + 2*sum(fx(2:n)) );    % midpoints    xm = (a+h/2:h:b);    % function values at midpoints    fx = feval( f, xm);    % composite Simpson rule    S(n) = S(n) + (2*h/3)*sum(fx);endfprintf('    n  &       T(n)        &        S(n)       &       int \n' )for n = 1:16  fprintf('  %3d  &  %14.13f  &  %14.13f  & %14.13f \n', n+1, T(n), S(n), int )end