scialb

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

function y=g(x)
    y=cos(x)
endfunction

n=100 //nb rectangles
clf
x=0:0.01:1
y=feval(x,f)
plot(x,y,'-r')
xgrid(3)

function animRect(f,a,b,n)
    disp(rect(f,a,b,n))
    ak=a:(b-a)/n:b
    yk=feval(ak,f)
    for i=1:n
        plot([ak(i),ak(i),ak(i+1),ak(i+1)],[0,yk(i),yk(i),0],'-m')
    end
endfunction

animRect(f,0,1,n)