function [T,Interval_Partition]=Romberg_Integral_Algorithm(fun,interval,num)

if  num==0;

   T=trapz(interval,feval(fun,interval));Interval_Partition=interval;

else

   [T_temp,I_temp]=Romberg_Integral_Algorithm(fun,interval,num-1);

   T=zeros(num+1,num+1);

   Interval_Partition=zeros(1,2*length(I_temp)-1);

   Interval_Partition(1:2:end)=I_temp;

   Interval_Partition(2:2:end-1)=(I_temp(1:end-1)+I_temp(2:end))/2;

   T(1:end-1,1:end-1)=T_temp;

   T(end,1)=trapz(Interval_Partition,feval(fun,Interval_Partition));

   for ind=2:num+1;

       T(end,ind)=4^(ind-1)/(4^(ind-1)-1)*T(end,ind-1)-1/(4^(ind-1)-1)*T(end-1,ind-1);

   end

end