close all

%% Analytical solutions a=1, b=2
format compact
format long
Io = [  2-sqrt(2)
        2^(3/2)-2
        (2^(5/2))/3 - 2/3
        (2^(7/2)-2)/5
    ]

ksiVec = [-3/2 -1/2 1/2 3/2];
j=1;
for ksi = ksiVec
    f=@(x) x.^ksi;

    % Bounds
    a=1;
    b=2;

    n = [10 20 40 80];  % Number of intervals
    h = (b-a)./n;       % Interval length.

    i = 1;
    for n = n
        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(j));
        i = i+1;
    end

    subplot(2,2,j)
    plot(h, Imidp, '*-b')
    title(['$\xi= $' num2str(ksi)], 'interpreter','latex')
    ylabel('I');
    xlabel('h');
    
    %%

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

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

    xVec = A \ bVec;

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

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

    legend('I_{mid}', 'I_{LS}')

    j=j+1;
end

