simp2.m

function q = simp2(y,h);
% function simp2(y,h);     Simpson's 1/3 quadrature rule - vector version
%
%   input: y - function vector to be integrated
%          h - interval width
%  output: integral in the given region

n = length(y);             % number of data points
if (n/2) == floor(n/2)     % check to see if even (don't want this!)
  fprintf('\nError: must have odd number of data points\n');
  break;
end

q = h/3 * (y(1) + 4 * sum(y(2:2:n-1)) + 2 * sum(y(3:2:n-2)) + y(n));