function f(x)
    return 1/(1+x^2)
end

ye=zeros(7);
he=zeros(7);
using Printf
for j=1:6
    N=2^j
    n=3
    a=0
    b=1
    h=(b-a)/(N*n)
    he[j]=h
    xi=a:h:b
    yi=f.(xi)

    w=zeros(N*n+1)
    for k=1:N
        w[k*n-2:k*n+1]+=[1;3;3;1]
    end

    s=w'*yi*3*h/8
    ye[j]=abs(s-atan(1))
    @printf("s=%.15e e=%.15e\n",s,s-atan(1))
end