function [ ] = romberg(aa, bb, kk, iter)
% Romberg Integral Method
% Coded by Nan.Xiao 2010-05-23
% Step 1. Calc TrapeStep
% Step 2. Calc Romberg Method
% Step 3. Extract Diagnal & Output
function [ ] = trapestep(a, b, k)
    format long
%     function y = f(x)
%         y=(2/sqrt(pi))*exp(-x);
%     end
%     Examples Can Also Be:
%     function y = f(x)
%         y=x*sin(x);
%     end
%     Or
    function y = f(x)
        y=x*sqrt(1+x.^2);
    end
k1=2^(k-1);
temp=zeros(k,k1);
mysum=zeros(k,1);
t=zeros(k+1,1);
for i=1:k
    for j=1:(2^(i-1))
    temp(i,j)=f(a+(((2*j-1)*(b-a))/(2^i)));
    end
end
for p=1:k
    mysum(p,1)=sum(temp(p,:));
end
t(1,1)=0.5*(f(a)+f(b));
for q=2:k+1
    t(q,1)=(0.5*t(q-1,1))+(((b-a)*mysum(q-1,1))/(2^(q-1)));
end
end
    format long
trapestep(aa, bb, kk);
diagonal=zeros(kk,1);
tt=zeros(kk,kk);
for v=1:kk
    tt(v,1)=t(v,1);
end
for r=2:kk
for s=2:r
    tt(r,s)=((4^(s-1)*tt(r,s-1))-tt(r-1,s-1))/(4^(s-1)-1);
end
end
save('tt.txt', 'tt', '-ascii', '-double'); % Save the Table
% Extract Table Diagonals
for w=1:kk
    diagonal(w,1)=tt(w,w);
end
disp('Romberg Table Diagonals:')
disp(diagonal);
% Output Final Result
cc=2;
while abs(diagonal(cc,1)-diagonal(cc-1,1))>=iter || abs(diagonal(cc+1,1)-diagonal(cc,1))>=iter
    cc=cc+1;
end
disp('The Final Result is:');
disp(diagonal(cc,1));
end