function int = sinx2(a,b,err)

%% uses simpsons rule to compute the approximate definite integral of
%% sin(x^2).
%% a and b are the enpoints of the interval on which the estimate is made
%% err is the tolerance. Stepsize is computed automatically.
%%initial estimate of fourth deriviative

srate=1e-2;
tint=a:srate:b-srate;                                                       % for derivative est
u=sin(tint.^2);
max4diff=abs(max(diff(diff(diff(diff(u./srate)./srate)./srate)./srate)))       % low res 4th derivative estimate

h=min([0.1,(err*180/((b-a)*max4diff))^(.25)]);   
N=uint32(floor((b-a)/h));                               % number of terms to compute
disp(['Number of rectangles used ',num2str(N)]);
t=a:h/2:b;                                      % should be 2N+1 entries
f=sin(t.^2);                                    % integrand
teven=t(1:2:length(t));                         % even (endpoint) samples (N+1)
todd=t(2:2:length(t)-1);                        % odd (midpoint) samples (N)

s=zeros(N,1);
for k=1:N
    s(k)=f(2*k-1)+4*f(2*k)+f(2*k+1);            % k-th rectangle sum
end
sum=0;
S=zeros(N,1);                                   % numerical integral
for l=1:N                                       % add next rectangle sum
    S(l)=sum+s(l);
    sum=S(l);
end
size(S)
S=h*S./3;                                       % to get right magnitude
int=S;                                          % output


plot(todd,sin(todd.^2),'r');
hold on;
plot(todd,int);
title(['Integrand (red) and Simpson approximation, N= ',num2str(N),' subintervals']);