
% This function produces the Romberg integration of the function "f" from 
% a to b. The output will produce a number that is within 10^-8 of the 
% exact answer given. After 12 complete iterations the integration is 
% stopped. 
function output = Romberg(f,a,b,exact)
f_b = feval(f,b);
f_a = feval(f,a);

R(1,1) = 1/2*(b-a)*(f_b + f_a);
approx = R;
i = 1;
while (abs(exact - approx) > 10^(-8))
    n = 2^i;
    h = (b-a)/n;
    R(i+1,1) = 1/2*R(i,1);
    sum = 0;
    for k = 1:2^(i-1)
        sum = sum + 1/(a + (2*k-1)*h);
    end
    R(i+1,1) = R(i+1,1) + h*sum;
    approx = R(i+1,1);
    for j = 1:i
        R(i+1,j+1) = R(i+1,j) + 1/(4^(j)-1)*(R(i+1,j)-R(i,j));
        approx = R(i+1,j+1);
    end
    i = i+1;
    if (i>12)
        fprintf('too many iterations')
        exit
    end
end

output = R


%function output = Trapezoid(a,b,f)
% % this function will produce an approximation to the function f from
% % a to b. The given inputs are the left boundry a and the right boundry
% % b and f is the vector of the function evaluated at uniformly spaced 
% % points. The integration is approximated by the composite trapezoidal 
% % rule.
% 
% n = length(f)
% h = (b-a)/n;
% if (n == 2)
%     output = h/2*(f(2)-f(1));
% else
%     output = h*(sum(f(2:length-1))
%     