-- Adaptive Simpson Integration
-- simpson(a, b, t)
--
-- integrates f(x) between a and b with error of t
--
-- Simpson algoritm that recursivly call itself on partition that need to be
-- divided into more partition to reduce the error.
--
-- by Jesse Costales rcostale@nmsu.edu
-- algorithm adopted from 'Numerical Mathematics and Computing',
-- 3rd ed. by W. Cheney & D. Kincaid.
-- (use Radians)

e(a,b) = (f(a) + 4*(f(h(a,b))) + f(b))*((b-a)/6)
g(a,b) = e(a,h(a,b)) + e(h(a,b),b)
h(a,b) = (a + b)/2

simpson(a,b,t) = g(a,b) when abs(e(a,b)-g(a,b)) <= 15*t,
                 simpson(a,h(a,b),t/2) + simpson(h(a,b),b,t/2) otherwise

-- define function
f(x) = sin(x)

s(x) = simpson(0,x,5e-3)
i(x) = -cos(x) + 1

Xmin =0; Xmax = 2*pi 
plot f(X)
plot s(X)
plot i(X)