kill(n);
decimales:7$ prec:10^decimales$
f(x):=exp(-x^2)$
a:0$ b:1$ m:5$ n:2*m$ h:(b-a)/n$
integral:integrate(f(x),x,a,b)$
x(i):=a+i*((b-a)/n)$
rectcentral(n):= makelist(rectangle([x(k-1),0],
         [x(k),f((x(k-1)+x(k))/2)]), k,1,n)$

sumariemann:0$
for i: 1 thru n do
    sumariemann:float(sumariemann+
                 ((b-a)/n)*f((x(i-1)+x(i))/2))$
trapezoidal:float((h/2)*(f(a)+f(b)))$
for i: 1 thru n-1 do
    trapezoidal:float(trapezoidal+h*f(x(i)))$
simpson:float((h/3)*(f(a)+f(b)))$
for i: 1 thru m do
simpson:float(simpson+4*h/3*f(x(2*i-1)))$
for i: 1 thru m-1 do
simpson:float(simpson+2*h/3*f(x(2*i)))$

print ("Riemann","   Trapecio"," Simpson")$
print (float(round(sumariemann*prec)/prec),
       float(round(trapezoidal*prec)/prec),
       float(round(simpson*prec)/prec))$
print ("El valor de la integral es ",
       float(round(integral*prec)/prec))$

load("draw")$
wxplot2d(f(x),[x,a,b],[legend,false])$
grafriemann(n):= apply(wxdraw2d, 
    append([xrange=[a,b], color=blue],
    rectcentral(n),[explicit(f(x),x,a,b)]) )$
grafriemann(n)$