function I=trapets(f,a,b,tol,I0)
m=(a+b)/2;
h=b-a;

fa=feval(f,a);
fm=feval(f,m);
fb=feval(f,b);

I1=h*(fa+fm)/4;
I2=h*(fm+fb)/4;
I=I1+I2;

if abs(I-I0)>tol;
   I1=trapets(f,a,m,0.5*tol,I1);
   I2=trapets(f,m,b,0.5*tol,I2);
   I=I1+I2;
else
    hold on
    plot([a a m m],[0 f(a) f(m) 0])
    plot([m b b],[f(m) f(b) 0])
    x=a:h/100:b;plot(x,f(x),'r')
    hold off
end
