simp1.m

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

n = length(y) - 1;         % number of regions
if (n/2) ~= floor(n/2)     % check to see if even
  fprintf('\nError: must have odd number of data points\n');
  break;
end

v  = 2 * ones(n+1,1);      % these commands form : 2 4 2 4 ... 4 2 4 2
v2 = 2 * ones(n/2,1);
v(2:2:n) = v(2:2:n) + v2;     

v(1) = 1;                  % this correct the end points
v(n+1) = 1;

q = y * v;                 % row vector * column vector = sum(y_i * v_i)
q = q * h/3;