
function demo_romberg

m=6;

N0=1;
N=N0*2.^(0:m-1);

example=1;

switch example
    case 1
 a=-1; b=1; f=@(x) exp(x); sol=exp(b)-exp(a);
    case 2
 a=-1; b=1; f=@(x) cos(x); sol=sin(b)-sin(a);
    case 3
a=0; b=pi/2; f=@(x) (x.^2+x+1).*cos(x); sol=-2+pi/2+pi^2/4;
end

% prima riga
for ii=1:m
    [x,w]=trapezi_composta(N(ii),a,b);
    fx=feval(f,x);
    I(ii,1)=w'*fx;
end


% righe successive
for k=1:m
    for ii=1:k
        I(ii,k+1)=NaN;
    end
    prev_col=I(k:m,k);
    curr_col=(2^(2*k)*prev_col(2:end)-prev_col(1:end-1))/(2^(2*k)-1);
    I(k+1:end,k+1)=curr_col;
end

format short e
errs=abs(I-sol);
h=(b-a)./N';
A=[h errs];
for ii=1:m
    for jj=1:m+1
        if isnan(A(ii,jj)) == 1
            fprintf('$$ &');
        else
        fprintf('$%1.1e$ &',A(ii,jj));
        end
    end
    fprintf(' \\\\\n \t');
end


fprintf('\n \n');








function [x,w]=trapezi_composta(N,a,b)

% FORMULA DEI TRAPEZI COMPOSTA.

% INPUT:
% N: NUMERO SUBINTERVALLI.
% a, b: ESTREMI DI INTEGRAZIONE.

% OUTPUT:
% x: NODI INTEGRAZIONE.
% w: PESI INTEGRAZIONE (INCLUDE IL PASSO!).

h=(b-a)/N;             % PASSO INTEGRAZIONE.
x=a:h:b; x=x';         % NODI INTEGRAZIONE.
w=ones(N+1,1);         % PESI INTEGRAZIONE.
w(1)=0.5; w(N+1)=0.5;
w=w*h;

