deff('y=f(x)','y=4/(1+x^2)');
t1=1/2*(f(0)+f(1));t=log(abs(t1-%pi));
for p=1:10;
n=2^(p-1);s=0;
for j=1:n
 s=s+f((2*j-1)/(2*n));
end
t1=1/2*(t1+1/n*s);t=[t,log(abs(t1-%pi))];
end
plot2d(0:10,t)
