In your solution the x-values are computed as follows:
h = (b-a)/n
x1 = a+1
x3 = x1 +2*h
x5 = x3 +2*h
...
This means rounding errors slowly accumulate.
It happens when (b-a)/n
is not representable as floating point.
If we instead compute xi
as a+ (i*(b-a))/n
you will get more accurate results.
This variant of your solution uses the above method to compute the xi
.
(define (simpson-integral3 f a b n)
(define h (/ (- b a) n))
(define (next i) (+ i 2))
(define (f* i) (f (+ a (/ (* i (- b a)) n))))
(* (/ h 3)
(+ (f a)
(* 4 (sum f* 1 next n))
(* 2 (sum f* 2 next (- n 1)))
(f b))))