function Z = srule(f,a0,b0,tol0)
% Z = srule(f,a0,b0,tol0)
% Subroutine for adaptive quadrature.
h = (b0 - a0)/2;
c0 = (a0 + b0)/2;
Fa = feval(f,a0);
Fc = feval(f,c0);
Fb = feval(f,b0);
S = h*(Fa + 4*Fc + Fb)/3;
S2 = S;
tol1 = tol0;
err = tol0;
Z = [a0 b0 S S2 err tol1];