After noticing the structure of the Simpson Integration Formula I realized that for interated functions this simplified somewhat. The two formulas are simpf1 , for Mandelbrot Set , and simpfn1 for a general function fn1 . Simpson integration is a low order approach , therefore most of the fractal structure of the iterated function appears to be preserved . simpf1(XAXIS) {; sciwise ; Simpson's integration applied to f1. c = Pixel s = 0 x0 = c x1 = (x0*x0) + c x2 = (x1*x1) + c : s = s + (x1 - x0)*( x2 + x1) x0 = x1 x1 = x2 x2 = (x2*x2) + c z = 0.5*abs(s) |z| < 40000 } simpfn1(XAXIS) {; sciwise ; Simpson's integration applied to f1. c = Pixel s = 0 x0 = c x1 = fn1(x0) x2 = fn1(x1) : s = s + (x1 - x0)*( x2 + x1) x0 = x1 x1 = x2 x2 = fn1(x2) z = 0.5*abs(s) |z| < 40000 }