
function demo_romberg2

n=6;

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

[r,h] = romberg(f,a,b,n);
ae=abs(r-sol);

A=[h' ae];

M=size(A,1);

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');


% Usage: r = romberg(f,a,b,n)
% Romberg integration
%
% Input:
% f - Matlab inline function 
% a,b - integration interval
% n - number of rows in Romberg tableau
%
% Output:
% r - Romberg tableau containing the computed values of the integral
%
% Examples:
% r=romberg(@sin,0,1,5);
% r=romberg(@myfunc,0,1,5);          here 'myfunc' is any user-defined function in M-file
% r=romberg(inline('sin(x)'),0,1,5);
% r=romberg(inline('sin(x)-cos(x)'),0,1,7);

function [r,h] = romberg(f,a,b,n)
h = (b - a) ./ (2.^(0:n-1));
r(1,1) = (b - a) * (f(a) + f(b)) / 2;
for j = 2:n
    subtotal = 0;
    for i = 1:2^(j-2)
        subtotal = subtotal + f(a + (2 * i - 1) * h(j));
    end
    r(j,1) = r(j-1,1) / 2 + h(j) * subtotal;
    for k = 2:j
        r(j,k) = (4^(k-1) * r(j,k-1) - r(j-1,k-1)) / (4^(k-1) - 1);
    end
end;