clc, clear

a = 0;
b = 1;
h = b-a;
T(1) = h/2*(a^(3/2)+b^(3/2));

%% 迭代得到 G0
for n = 1:10
    
    hn = h/2^(n-1);
    xn = a + ( 0:2^(n-1)-1 )*hn + hn/2;
    T(n+1) = T(n)/2 + hn/2 * sum( xn.^(3/2) );
    
end

%% Romberg 方法
G = zeros(size(T,2));
G(:,1) = T;  

for m = 2: size(G,2)
    
    G(m:end,m) = ( 4^(m-1)*G(m:end,m-1) -G(m-1:end-1,m-1) )/(4^(m-1)-1);
    if abs(G(m,m)-G(m-1,m-1)) < 1e-7
        fprintf('\n  G(%d) = %.15f\n |G(%d)-G(%d)|=%s\n',m,G(m,m),m,m-1,abs(G(m,m)-G(m-1,m-1)) );
        break;
    end
    
end