%?龍貝格一維/二維積分公式。積分值=rombeg（‘函數‘XminXmax[YminYmax][精度隱含精度1e-8]）
%?二維情況下可做變上下限積分（必須是匿名函數形式）.此程序積分區間和函數可為矩陣
%?但兩者不能同時為矩陣。變上下限積分時只能y變量是矩陣計算效果等同x有矩陣情況
function?I=rombeg（fXminXmaxYminYmaxtol）
if?~isa（f‘function_handle‘）&&~isa（f‘inline‘）
???f=inline（f）;
end
switch?nargin
????case?3?%一維積分
?????????I=jf1（fXminXmax）;?%龍貝格一維積分計算函數
????case?4%龍貝格一維積分計算函數
????????I=jf1（fXminXmaxYmin）;
????case?5??%二維積分
????h=Xmax-Xmin;
????F1=@（y）?f（Xminy）;
????F2=@（y）?f（Xmaxy）;
?????if?isa（Ymin‘function_handle‘）
???????c1=Ymin（Xmin）;
???????c2=Ymin（Xmax）;
?????else?c1=Ymin;c2=Ymin;
?????end
?????if?isa（Ymax‘function_handle‘）??
?????????d1=Ymax（Xmin）;
?????????d2=Ymax（Xmax）;
?????else?d1=Ymax;d2=Ymax;
?????end
????s1=jf1（F1c1d1）;
????s2=jf1（F2c2d2）;
????V（::1）=（s1+s2）.*h/2;
????n=1;
????k=1;
????while?k<=10
????????n=2*n;
????????h=h/2;
????????k=k+1;
????????sum=zeros（size（s1））;
????????for?p=1:n/2
????????????F=@（y）?f（Xmin+（2*p-1）*hy）;
?????????????if?isa（Ymin‘function_handle‘）
????????????c=Ymin（Xmin+（2*p-1）*h）;
?????????????else?c=Ymin;
?????????????end
?????????????if?isa（Ymax‘function_handle‘）??
????????????d=Ymax（Xmin+（2*p-1）*h）;
?????????????else?d=Ymax;
?????????????end
????????????s=jf1（Fcd）;
????????????sum=sum+s;
????????end
????????V（::k）=V（::k-1）/2+（sum.*h）;
????????m=k;
??????????while?1%龍貝格遞推
????????????V1（::m-1）=（4*V（::m）-V（::m-1））/3;
????????????m=m-1;
?????????????if?m==1
??????????????break
?????????????else?V2（::m-1）=（8*V1（::m）-V1（::m-1））/7;
??????????????m=m-1;
?????????????end
?????????????if?m==1
?????????????????break
???????????????else?V3（::m-1）=（16*V2（::m）-V2（::m-1））/15;
?????????????????????m=m-1;
?????????????????????break
?????????????end
??????????end
??????????????if?k==4&&（max（（max（abs（real（V3（::1）-V2（::2）））））））<=1e-8&&（max（（max（abs（imag（V3（::1）-V2（::2）））））））<=1e-8%精度控制
?????????????????I=V3（::m）;
????????????????break
?????????????elseif?k>4&&（max（max（abs（real（V3（::m）-V3（::m-1））））））<=1e-8&&（max（（max（abs（imag（V3（::1）-V2（::2）））））））<=1e-8
?????????????????I=V3（::m）;
??????????????break
??????????????end
????end
case?6
????h=Xmax-Xmin;
????F1=@（y）?f（Xminy）;
????F2=@（y）?f（Xmaxy）;
?????if?isa（Ymin‘function_handle‘）
???????c1=Ymin（Xmin）;
???????c2=Ymin（Xmax）;
?????else?c1=Ymin;c2=Ymin;
?????end
?????if?isa（Ymax‘function_handle‘）??
?????????d1=Ymax（Xmin）;
?????????d2=Ymax（Xmax）;
?????else?d1=Ymax;d2=Ymax;
?????end
????s1=jf1（F1c1d1）;
????s2=jf1（F2c2d2）;
????V（::1）=（s1+s2）.*h/2;
????n=1;
????k=1;
????while?k<=10
????????n=2*n;
????????h=h/2;
????????k=k+1;
?????