function [Q,R] = romberg(a,b,n)
  
h = b-a;
R(1,1) = (h/2)*(f(a)+f(b));

for i=1:n
  
  h = h/2;
  sum = midtpunkt(a,b,i);
  R(i+1,1) = R(i,1)/2+sum;
  for j=1:i
    R(i+1,j+1) = R(i+1,j)+(R(i+1,j)-R(i,j))/(4^(j)-1);
  end
  
end
Q = R(n+1,n+1);
R
  
% Funksjon som regner ut de indre punktene i rekursiv trapes (se s.197)
function U = midtpunkt(a,b,n)

h = (b-a)/2^n;
nin = 2^(n-1);
x = a+h*(2*(1:nin)-1);
y = f(x);
U = h*sum(y);
  
function y = f(x)

y = (1+x).^(-1);
%y = exp(x);
%y = (1+x.^2).^(-1);
%y = sqrt(x);