clear all 
close all
a=0;%end points of intergration region
b=1;
m=0;% value for m  (using the notation from the textbook)
hm=(b-a)/2^m;%from m we calculate h
f=@(x) sin(x); %integrand


format long%print lots of decimal places
RapxI=Roberg(f,a,b,m,3)%aproximate the integral

exactI=1-cos(1)% exact value for integral

%% error as a function of m (or equivilantly h)
mmin=0;
mmax=12;
mlist=[mmin:1:mmax];% list of ms
hlist=(b-a)./2.^mlist;%list of hs
error=0*mlist;
for i=1:length(mlist);
error(i)=1-cos(1)-Roberg(f,a,b,mlist(i),3);% list of errors
end

figure%plot error vs h on log scale
loglog(hlist,abs(error),'o--','LineWidth',2)
hold on
plot(10.^(-[0:0.1:1]),10.^(-8*[0:0.1:1]-10),'r-.','LineWidth',2)
plot([10^(-4),1],10^(-16)*[1,1],'r:','LineWidth',2)
xlabel('$h$','interpreter','latex','FontSize', 18)
ylabel('$|$error$|$','interpreter','latex', 'FontSize', 18)
print( '-depsc2', 'Q2err' )


function IntVal=Roberg(f,a,b,m,k)% Romberg integration function
% f function handle of the integrand
% a lower limit of integral
% b upper limit of integral
% m the inital number of divisions will be 2^m
% k the largest number of divisions will be 2^(m+k)
Ts=0*[1:1:k+1];% list of T's 

%% Calculate T_{m,0} using formula (2.1) in solutions
hm=(b-a)/2^m;%h(m)
Ts(1)=(sum(feval(f,[a:hm:b]))-(feval(f,a)+feval(f,b))/2)*hm;

%%
for j=1:k% calculating the T_{m+j,0} for j=1,...,k using eq (2.3) in solutions
hj=hm/2^j;%h(j+m)
Ts(j+1)=Ts(j)/2+hj*sum(feval(f,a+(2*[1:1:2^(m+j-1)]-1)*hj));%eq (2.3) 
%note Ts(j+1)=T_{m,j}
end


for j=1:k%Richardson extrapolation
    oldTs=Ts;% list oldTs=[T_{m-j+1,j-1}, T_{m+2-j,j-1},...,T_{m+k,j-1}]
    Ts=0*[1:1:length(oldTs)-1];Ts;
    % Ts will be [T_{m-j+2,j}, T_{m+3-j,j},...,T_{m+k,j}]
    for jj=1:k+1-j% could probably speed this up by vectorising this
        Ts(jj)=(2^(2*j)*oldTs(jj+1)-oldTs(jj))/(2^(2*j)-1);%equation (2.2)
    end
    
end

IntVal=Ts;
end
