function [I]=simp13(t,y)
% This function produces an integral estimate for data given
% by a pair of independent/dependent data row vectors.
% It is most appropriate to use this function with 3 or more data points, but
% values will be returned if there are at least 2. In that case the
% trapezoid rule is employed.
% The first and last interval integrals are calculated by using
% the best-fit quadratic (i.e. The Simpson 1/3 rule).
% Interior interval integrals are calculated by averaging
% the integrals of two quadratics.
% If the t entries are gigantic or very small MATLAB might complain (though
% it will give an answer) about a poorly conditioned matrix, since the
% backslash operation on a matrix composed of powers of these entries is
% employed to find the cubic coefficients.
% Try 10^n*simp13(t/10^n,y) for appropriate large (positive or negative) n
% to silence this complaint.
if length(t)~=length(y) error('the two data vectors must have equal length'); end
if unique(sort(t))~=t error('the domain data vector must have distinct entries listed in increasing order'); end
d=length(t); % Number of data points
n=d-1; % Number of data intervals
if d==2 I=(t(2)-t(1))*(y(1)+y(2))/2; end
if d==3
A=[t(1)^2 t(2)^2 t(3)^2 ; t(1) t(2) t(3) ; 1 1 1]'\y(1:3)';
B=[1/3*A(1) 1/2*A(2) A(3) 0];
I=polyval(B,t(3))-polyval(B,t(1));
end
if d>3
M=zeros(d-2,d-1);
for k=0:d-3;
A=[t(k+1)^2 t(k+2)^2 t(k+3)^2 ; t(k+1) t(k+2) t(k+3) ; 1 1 1 ]'\y(k+1:k+3)';
B=[1/3*A(1) 1/2*A(2) A(3) 0];
M(k+1,k+1)=polyval(B,t(k+2))-polyval(B,t(k+1));
M(k+1,k+2)=polyval(B,t(k+3))-polyval(B,t(k+2));
end
for j=1:d-1
deltaI(j)=sum(M(:,j))/2;
end
deltaI(1)=2*deltaI(1);
deltaI(d-1)=2*deltaI(d-1);
I=sum(deltaI(1:d-1));
end