function int=intspline(x,y)
%
%  function generates a cubic spline and integrates it
%  from x(1) to x(n) using Simpson's rule
%
n=size(x,1);
[a,b,c,d]=genspline(x,y,n);
int=0;
for i=1:n-1
   f1=evalspline(a,b,c,d,x(i),x,n);
   f2=evalspline(a,b,c,d,.5*(x(i)+x(i+1)),x,n);
   f3=evalspline(a,b,c,d,x(i+1),x,n);
   int=int+(x(i+1)-x(i))/6*(f1+4*f2+f3);
end
int
