```
(bind-func nintegrate
"integrate using Simpsons rule - N must be odd"
(lambda (f:[double,double]* xmin xmax N)
(let ((sumeven:double 0.0)
(sumodd:double 0.0)
(x 0.0) (n 0)
(interval (/ (- xmax xmin) (i64tod (- N 1)))))
(set! n 2)
(while (< n N)
(set! x (+ xmin (* interval (i64tod (- n 1)))))
(set! sumodd (+ sumodd (f x)))
(set! n (+ n 2)))
(set! n 3)
(while (< n N)
(set! x (+ xmin (* interval (i64tod (- n 1)))))
(set! sumeven (+ sumeven (f x)))
(set! n (+ n 2)))
(* (+ (* sumeven 2.0) (* sumodd 4.0)
(f xmin) (f xmax))
(/ interval 3.0)))))
```

Back to Index