% -------------------------------------------------------------------------
% ----------------------------龍貝格積分算法--------------------------------
% 函數  In=Romberg(func,a,b)
% 輸入參數: func   被積函數的名稱,如'sin' 
%          a,b    積分區間
% 輸出結果: In     積分結果,在命令窗口顯示ans=____
%    >> I = Romberg( 'sin', 0, pi/2)
% 編程思路:1.設計一個矩陣rb中,根據算法得出矩陣元素之間的關系
%         2.積分限由(a,b)--->(0,1)  便于尋找規律進行編程
%         3.也有可能會出現計算很淺（行數），就終止的情況，這會出現很大誤差！?。?%           譬如：f(0.5)=0.5*[f(0)+f(1)],第二行就停了。
%           為避免此種情況，讓程序先自動運行jingdu行，再判斷誤差，用大計算量換取可靠。
% -------------------------------------------------------------------------
function In=Romberg(func,a,b)
wucha=1e-8;
jingdu=6;
rb(1,1)=(b-a)*(feval(func,a)+feval(func,b))/2;  
i=1;
while i<=jingdu||abs(rb(i,i)-rb(i,i-1))>wucha  %第一步判斷條件時i＝1，但是執行不到abs(rb(i,i)-rb(i,i-1))>wucha!!!
    i=i+1;
    temp=0;
    for k=1:2:2^(i-1)
        temp=temp+(b-a)*feval(func,a+(b-a)*k/(2^(i-1)));
    end
    rb(i,1)=rb(i-1,1)/2+temp/(2^(i-1));%T
    for j=2:i
        rb(i,j)=(-rb(i-1,j-1)+(4^(j-1))*rb(i,j-1))/(4^(j-1)-1);%S T C R
    end
end
In=rb(i,i);