%% Example 3.2
%  Romberg integration


format long;


%% I_1, I_2, I_3, I_4
I = zeros(4,1);
h = zeros(4,1);

for i = [1:4]
    N    = 2^i;
    h(i) = (pi-1)/N;
    X    = [1:h(i):pi];
    f    = sin(X)./(2*X.^3);
    
    for j = [2:N]
        I(i) = I(i) + f(j);
    end;
    I(i) = I(i) + 1/2 * (f(N+1) + f(1));
    I(i) = I(i)*h(i);
end;


%% I_12, I_23, I_24

I_12 = (4*I(2)-I(1))/3;
I_23 = (4*I(3)-I(2))/3;
I_34 = (4*I(4)-I(3))/3;

%% I_123, I_234

I_123 = (16*I_23-I_12)/15;
I_234 = (16*I_34-I_23)/15;


%% I_1234

I_1234 = (64*I_234-I_123)/63








