h:= ('b'-'a')/'n';
X:= i -> a + i*h;
between:= (x,L,R) -> 0 < evalf(x - L) and 0 < evalf(R - x);
upper:= i -> max(f(X(i - 1)), f(X(i)),
     op(map(f, select(between, cp, X(i - 1), X(i)))));
lower:= i -> min(f(X(i - 1)), f(X(i)),
     op(map(f, select(between, cp, X(i - 1), X(i)))));
rsum   := g -> h * sum('g(i)', i=1..n);
rect   := (x1,x2,y) -> [[x1, 0], [x1, y], [x2, y], [x2, 0]];
rbox:= proc(g) local ii;
       plot({ f(x), 
             [ seq(op(rect(X(ii-1),X(ii),g(ii))),ii=1..n)]},
            x = a..b, args[2..nargs]);
       end;