http://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method
double simpson_rec( double a, double b, double S, double fa, double fb, double fc, int bottom, double epsilon){ double c = (a+b)/2, h = b-a; double d = (a+c)/2, e = (c+b)/2; double fd = f(d), fe = f(e); double SL = (h/12)*(fa + 4*fd+fc); double SR = (h/12)*(fc + 4*fe+fb); double S2 = SL + SR; if( bottom <= 0 || fabs(S2 - S ) <= 15*epsilon){ return S2 + (S2-S)/15; } return simpson_rec(a,c,SL,fa,fc,fd, bottom-1, epsilon/2) + simpson_rec(c,b,SR,fc,fb,fe, bottom-1, epsilon/2); } double simpson( double a, double b, double epsilon, int rec){ double c = (a+b)/2, h = b-a; double fa = f(a); double fb = f(b); double fc = f(c); double S = (h/6) * (fa+4*fc+fb); return simpson_rec(a,b,S,fa,fb,fc,rec,epsilon); }