Below are various implementations of an adaptive integration method (based on the trapezoid rule) using Mathematica's recursive definition feature.  The Timing[] function is used to compare the computational efficiency of each implementation.

Set up

Clear[f, t1, t2, adaptive, adaptive1, adaptlist]

f[x_] = Exp[-x^2]

exact = NIntegrate[Exp[-x^2], {x, 0, 2}]

Version 0: ignore function evaluation efficiency

t1[a_, b_] := (b - a)/2*(f[a] + f[b]); 
    t2[a_, b_] := (b - a)/2/
          2*(f[a] + 2  f[(a + b)/2] + f[b]);

adaptive[a_, b_, 
          tol_] := 
If[
Abs[
              t1[a, b] - t2[a, b]] < 3  tol, 
t2[a, 
            b], 
adaptive[a, (a + b)/2, tol/2] + 
            adaptive[(a + b)/2, b, 
              tol/2]
];

approx = adaptive[0. , 2. , 0.01]

exact - approx

Version 0.1: list version

This version builds a list of lists.  Each inner list  is an ordered triple that gives the endpoints of a terminal interval in the recursion and the contribution to the approximation from that interval.  We can sum the contributions to get the total approximation.  We can use the endpoints to make a plot showing where the function is evaluated.

Trapezoid

adaptivelist[a_, b_, tol_] := 
If[
          Abs[t1[a, b] - t2[a, b]] < 3  tol, 
          approxlist = 
            Append[approxlist, {a, b, 
                t2[a, b]}], 
adaptivelist[a, (a + b)/2, 
            tol/2]; adaptivelist[(a + b)/2, b, tol/2]];