% Example 11.7
% Romberg's Method
clearvars
clc

f = @(x) 2 + 2.*x + x.^2 + sin(2.*pi().*x) + cos(2.*pi().*x./0.5);
a = 0;
b = 1.5;
n1 = 1;

% Calculate the exact integral
Iexact = integral(f,a,b);

% TRAPEZOIDAL RULE
disp(' ')
disp('TRAPEZOIDAL RULE')
disp('--------------------------------------------------')

trapmat = zeros(8,5);
for i=1:8
    n = 2^(i-1)*n1;
    h = (b-a)/n;
    IT = trap_func(f,a,b,n);
    et = abs(Iexact - IT);
    trapmat(i,:) = [n h IT Iexact et];
end

fprintf('\n  n          h       IT        Itrue           |E|\n')
fprintf('%3g %11.7g %10.6f %10.6g %14.7f\n',transpose(trapmat));



% ROMBERG's METHOD
disp(' ')
disp(' ')
disp('ROMBERG METHOD')
disp('--------------------------------------------------')

% Initialize values
maxit = 5;  %Stop at depth k = 5

count = 0;
er = 1;

IR = zeros(1);
k = 1;

while (count<maxit)
    count = count + 1;
    for i = 1:count
        IR(i,1) = trap_func(f,a,b,2^(i-1)*n1);
    end
    if(count>1)
        for k = 2:count
            for j = 1:count-k+1
                IR(j,k) = (4^(k-1)*IR(j+1,k-1) - IR(j,k-1))/(4^(k-1)-1);
            end
        end
    end
    er = Iexact - IR(1,k);
   
    outmat(count,1:2) = [count, abs(er)];    
end

Ifinal = IR(1,k);

 disp(['Romberg Table for k = ',num2str(count)])
 disp(' ')
 disp(IR)    


disp(['The estimate of the integral is ',num2str(Ifinal),' after ',num2str(count),' iterations'])
disp(' ')
disp(['The error estimate corresponding to each iteration is:'])
fprintf('\n  #          |E|\n')
fprintf('%3g %14.7f\n',transpose(outmat));