clear all;

a = 0;
b = pi/2;

Iex = exp(b.*sin(4*b))-exp(a.*sin(4*a))+b.^3-a.^3;

NMAX  = 10;

for N=2:2:NMAX*2
h=(b-a)/N;
x=linspace(a,b,N+1);
Neval = N+1;
y=funz(x);
I=h*(y(1)+2*sum(y(2:N))+y(N+1))/2;
ErrTrap(Neval) = abs(I-Iex);
end

for N=1:NMAX
h=(b-a)/N;
x=linspace(a,b,2*N+1);
Neval = 2*N+1;
y=funz(x);
I=h*(y(1)+2*sum(y(3:2:2*N-1))+4*sum(y(2:2:2*N))+y(2*N+1))/6;
ErrSimp(Neval) = abs(I-Iex);
end

nodi

for N=3:2:NMAX*2+1
Neval = N;
x = a+(b-a)/2*(GaussLegendre(N).x+1);
y=funz(x);
I = (b-a)/2*GaussLegendre(N).w*y';
ErrGL(Neval) = abs(I-Iex);
end

file = fopen('tabella1.txt','w');
for N=3:2:NMAX*2+1
fprintf(file,'%d\t%.4E\t%.4E\t%.4E\n',N,ErrTrap(N),ErrSimp(N),ErrGL(N));
end
fclose(file);