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

N=5
n=3
a=0
b=1
h=(b-a)/(N*n)

xi=a:h:b
yi=f.(xi)

using Plots
plot(xi,yi,marker=2)

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