close all
clear

f=@(x)sin(x.^2)
n=[17:20]
a=0; b=100;

syms t
F = sin(t.^2)
Fint = int(F,[0 100])
Io = vpa(Fint)

subints=2.^n
h = (b-a)./subints;

i = 1
for n = subints
    x = linspace(a, b, n+1);
    Imidp(i) = sum(f((x(1:n)+x(2:n+1))./2).*h(i));
    err(i) = abs(Imidp(i) - Io);
    i=i+1;
end

    alpha = log ((Imidp(2) - Imidp(3)) / (Imidp(3) - Imidp(4)));
    alpha = alpha / ( log ( h(2) / h(3) ) );

    A = [ ones(length(h)-1,1) (h(2:end).^alpha)' ];
    bVec = Imidp(2:end)'; % skip h0

    xVec = A \ bVec;

    hold on
    hcont = linspace(h(length(h)-1)/100, h(1), 100);
    plot(hcont, xVec(1) + xVec(2).*hcont.^alpha, '-r')
    grid on

    Io_LS = xVec(1)
    abs(xVec(1) - Io)





plot(h, Imidp, '*-')
ylabel('err(h)')
xlabel('h')
grid on
title('$f(x) = \sin(x^2)$', 'interpreter','latex')
%Io = sqrt(2*pi)/4; % Analytical solution for a=0, b=+infty