function [R, M]=romberg(a,b,tol,maxit)
%[R, M]=Romberg(a,b,tol,maxit)
%
%Romberg integration uses trap1step and myf
global count;

k=1;
Ro=[(b-a)*(myf(a)+myf(b))/2];   %Trap on 1 subinterval
%count=count+2;
M=Ro;
done=0;
while ((done ~= 1)&(k < maxit))
    k=k+1;
    Rn=[trap1step(a,b,2^(k-1),Ro(1))];
    for j=1:k-1
        Rn=[Rn Rn(j)+(Rn(j)-Ro(j))/(4^j-1)];
    end
    if (abs(Rn(k)-Ro(k-1)) < tol)
        done=1;
    end
    Ro=Rn;
    M=[[M zeros(k-1,1)]; Ro];
end
if (done ~= 1)
    msg='maximum number of itterations reached'
end
R=Ro(k);