Simpson's Rule

/* User input comes first */
f: x^2*sin(x)  $/* Function f(x) */
a: 1   $/* Interval upper bound */
b: 3   $/* Interval lower bound */
N: 100   $/* Number of subintervals (must be even) */
/* End of user input */

if ( mod(N,2) # 0 ) then ( N: N+1 )$
h: float( b - a ) / N$
sum: ( ev( f, numer, x=a ) + ev( f, numer, x=b ) )$
for i: 1 thru N-1 step 2 do (
   sum: sum + 4 * ev( f, numer, x=a+i*h )
)$
for i: 2 thru N-2 step 2 do (
   sum: sum + 2 * ev( f, numer, x=a+i*h )
)$
print( "Simpsons rule approximation is: ", h*sum / 3.0 )$