function S=rect(f,a,b,n)
    S=0
    ai=a
    for i=1:n
        S=S+(b-a)/n*f(ai)
        ai=ai+(b-a)/n
    end
endfunction

function y =f(x)
    y=x**3
endfunction

n=400
clf
x=0:0.001:1
ak=0:(1/n):1
y=feval(x,f)
yk=feval(ak,f)
plot(x,y,'-r')

for i=1:n
    plot([ak(i),ak(i),ak(i+1),ak(i+1)],[0,f(ak(i)),f(ak(i)),0],'-b')
end
xgrid(3)

disp(rect(f,0,1,n))