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

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

function animTrap(f,a,b,n)
disp(trap(f,a,b,n))

ak=a:((b-a)/n):b
yk=feval(ak,f)

for i=1:n
    if i<n/3 then
    plot([ak(i),ak(i),ak(i+1),ak(i+1)],[0,yk(i),yk(i+1),0],'-b')
    elseif i<n-n/3
    plot([ak(i),ak(i),ak(i+1),ak(i+1)],[0,yk(i),yk(i+1),0],'-w')
    else
    plot([ak(i),ak(i),ak(i+1),ak(i+1)],[0,yk(i),yk(i+1),0],'-r')
    end
end
endfunction

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

animTrap(f,0,1,n)