function [t2,M]=contra(f,a,b,tol)

%Input      - f is the integrand 
%              - a and b are upper and lower limits of integration
%              - tol is the tolerance
%Output   - t2 is the trapezoidal rule sum
%              - M is the number of intervals

%If f is defined as an M-file function use the @ notation
% call t2=contra(@f,a,b,tol).
%If f is defined as an anonymous function use the
% call t2=contra(f,a,b,tol).

%  NUMERICAL METHODS: Matlab Programs
% (c) 2007 by Xie Liling
%  Complementary Software to accompany the textbook:
%  Information and Computing Science: A Laboratory Course

h=b-a; M=2;
t1=(feval(f,a)+feval(f,b))*h/2;
h=h/2;
x=a+h;
t2=t1/2+feval(f,x)*h;

while (abs(t2-t1)>tol)
    M=2*M;
    h=h/2;
    t1=t2;
    s=0;
    for k=1:M/2
        x=a+h*(2*k-1);
        s=s+feval(f,x);
    end
    t2=t1/2+h*s;
end