function integ=SimpsonsRule(f,a,b)

%integ=SimpsonsRule(f,a,b)
%f is a row vector of function values.
%The independent variable must take equally spaced values
%on the interval [a,b], an odd number of values in all.
%integ is the integral of f from a to b by Simpson's Rule.
%f can also be a rectangular matrix, in which case each
%row contains function values, and the integral of each
%row from a to b is returned:  integ is then a column of values.

sizef=size(f); %size() returns [rows,columns]
rowlength=sizef(2); %number of columns
if(~mod(rowlength,2)) %if that is not odd (i.e., if it is even)
    'Integral not computed. Must supply an ODD number of values' 
else
    wt=2*mod(1:rowlength,2)+4*(~mod(1:rowlength,2)); %2 in odd places, 4 in evenplaces   
    wt(1)=1; %1 at the left end
    wt(rowlength)=1;  % and the right end
    wt=(b-a)/(rowlength-1)*wt/3; %wt is a row vector:
    integ=f*wt'; %the matrix f multiplies the column vector wt' 
end