function [I]=cumsimp38(t,y)

% This function produces the cumulative integral estimate for data given 
% by a pair of independent/dependent data row vectors.
% The output is a vector of cumulative integrals up to the data domain values.

% It is most appropriate to use this function with 4 or more data points, but 
% values will be returned if there are at least 2. Consider using trapz if you
% have 2 or 3 data points.

% The first and last interval integrals are calculated by using 
% the best-fit cubic (i.e. The Simpson 3/8 rule).
% The second and penultimate interval integrals are calculated by averaging 
% the integrals of two cubics
% Interior interval integrals are calculated by averaging 
% the integrals three cubics.

% 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*cumsimp38(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(1)=0; I(2)=(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(1)=0;
I(2)=polyval(B,t(2))-polyval(B,t(1));
I(3)=polyval(B,t(3))-polyval(B,t(1));
end

if d==4
A=[t(1)^3 t(2)^3 t(3)^3 t(4)^3 ; t(1)^2 t(2)^2 t(3)^2 t(4)^2 ; t(1) t(2) t(3) t(4) ; 1 1 1 1]'\y';
B=[1/4*A(1) 1/3*A(2) 1/2*A(3) A(4) 0];
I(1)=0;
I(2)=polyval(B,t(2))-polyval(B,t(1));
I(3)=polyval(B,t(3))-polyval(B,t(1));
I(4)=polyval(B,t(4))-polyval(B,t(1));
end

if d>4

M=zeros(d-3,d-1);

for k=0:d-4;

A=[t(k+1)^3 t(k+2)^3 t(k+3)^3 t(k+4)^3 ; t(k+1)^2 t(k+2)^2 t(k+3)^2 t(k+4)^2 ; t(k+1) t(k+2) t(k+3) t(k+4) ; 1 1 1 1]'\y(k+1:k+4)';
B=[1/4*A(1) 1/3*A(2) 1/2*A(3) A(4) 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));
M(k+1,k+3)=polyval(B,t(k+4))-polyval(B,t(k+3));
end


for j=1:d-1
deltaI(j)=sum(M(:,j))/3;
end

deltaI(1)=3*deltaI(1);
deltaI(2)=3*deltaI(2)/2;
deltaI(d-1)=3*deltaI(d-1);
deltaI(d-2)=3*deltaI(d-2)/2;

for j=1:d-1
I(j)=sum(deltaI(1:j));
end
I=[0, I];


end
