

close all % Close all open figure windows.
% let ksi = -3/2; -1/2; 1/2; 3/2
% ksi = [-3/2 -1/2 1/2 3/2];  
ksi = [3/2];
f=@(x)x.^ksi; % The function we evaluate.

a=1; b=2; % integral upper and lower bounds
n=10;
h=(b-a)/n;

h0 = (b-a)/n;
h1 = h0/2
h2 = h0/4
h3 = h0/8

hvec = [h0 h1 h2 h3]

for evalConVal = ksi
    i = 1;
    f = @(x) x.^evalConVal;
    controlValue = integral(f, a, b); % Accurate control integral for error computation.
    I0(i) = controlValue
    i = i+1;
end   
%% Middle Riemann sum (midpoint rule)
figure
for currentKsi = ksi
    i = 1;
    for n=[10 20 40 80]
        x = linspace(a,b,n+1);
        h = (b-a)/n;
        f = @(x) x.^currentKsi;
        controlValue = integral(f, a, b); % Accurate control integral for error computation.
        Imidp(i) = sum(f((x(1:n)+x(2:n+1))./2).*h)
        err(i) = abs(Imidp(i) - controlValue);
       
        i = i+1;
    end
    err;
    plot(fliplr(hvec), fliplr(err))
    hold on
    grid on
end
xlabel('h');ylabel('err(h)');

title(['Midpoint Rule. n= ', num2str(n), ' err= ', num2str(err)]);




%% Solving linear system of equations 
%  to obtain order of convergence alpha

alpha = log ((Imidp(1) - Imidp(2)) / (Imidp(2) - Imidp(3)));
alpha = alpha / ( log ( h1 / h2 ) )

A = [
    1 h0^alpha;
    1 h1^alpha;
    1 h2^alpha;
    1 h3^alpha;
    ] 
bVec = Imidp(1:end)' % skip h0

xVec = A \ bVec


%% plot least squares curve
hold off
figure
plot(fliplr(hvec), fliplr(Imidp))
hold on


hcont = linspace(h3, h0, 100);
LScurve = xVec(1)+xVec(2).*hcont.^alpha;
plot(hcont, LScurve)

figure
haxis = linspace(10e-3, 1, 100);
e = exp(log(xVec(2)) + alpha*log(haxis));
plot(haxis, e)