simp0.m

function q = simp0(y,h);
% function simp0(y,h);     Simpson's 1/3 quadrature rule - non-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

summat = 0;
for i = 2:2:n-1
  summat = summat + 4 * y(i);
end
for i = 3:2:n-2
  summat = summat + 2 * y(i);
end
summat = summat + y(1) + y(n);

q = summat * h/3;