--  Simpson integration

-- simpson(x1, x2, steps)  integrates f(x) between x1 and x2
-- by Matthieu Contensou   mc@astro.caltech.edu

g(x1, x2, n, s) = f(x1 + (x2-x1) * (n-1) / s) 

simpson(x1, x2, steps) = (sum(g(x1, x2, n, steps), n, 4, steps-2) +
  3/8*(g(x1, x2, 1, steps)   + g(x1, x2, steps+1, steps)) +
  7/6*(g(x1, x2, 2, steps)   + g(x1, x2, steps, steps)) +
  23/24*(g(x1, x2, 3, steps) + g(x1, x2, steps-1, steps))) *
  (x2 - x1) / steps



f(x) = x^2
simpson(0, 1, 500):0.333