function y = open1dK(f,a,b,m,n,r)

% This is a set up for 1d open romberg integration 
% f is the function
% a and b are the limits of integration
% n is the number of partitions in each subinterval
% m is the number of subintervals
% r is used to calculate the length of each subinterval
r1(1) = b;
for i = 1:m
    %r1(i+1) = 1/r^i;
    %r1(i+1) = 1/(i*r);
    %r1(i+1) = 1/sqrt(i*r);
    r1(i+1) = 1/(i*r)^(1/3);
end

for i = 1:m
    for k = 1:n
        %ank1(i,k) = (r1(i) - r1(i+1))/n;
        ank1(i,k) = 2*k*(r1(i) - r1(i+1))/(n^2 + n);
        %ank1(i,k) = 6*k^2*(r1(i)-r1(i+1))/((n*(n+1)*(2*n+1)));
        %ank1(i,k) = 4*k^3*(r1(i)-r1(i+1))/(n^2*(n+1)^2);
    end
end
for i = 1:m
    for k = 1:2*n
        %ank2(i,k) = (r1(i) - r1(i+1))/n;
        ank2(i,k) = k*(r1(i) - r1(i+1))/(2*n^2 + n);
        %ank2(i,k) = 6*k^2*(r1(i)-r1(i+1))/((2*n*(2*n+1)*(4*n+1)));
        %ank2(i,k) = 4*k^3*(r1(i)-r1(i+1))/(4*n^2*(2*n+1)^2);
    end
end
for i = 1:m
    for k = 1:4*n
        %ank3(i,k) = (r1(i) - r1(i+1))/n;
        ank3(i,k) = k*(r1(i) - r1(i+1))/(8*n^2+2*n);
        %ank3(i,k) = 6*k^2*(r1(i)-r1(i+1))/((4*n*(4*n+1)*(8*n+1)));
        %ank3(i,k) = 4*k^3*(r1(i)-r1(i+1))/(16*n^2*(4*n+1)^2);
    end
end

for i = 1:m
    right1(i,1) = r1(i);
    left1(i,n) = r1(i+1);
    for k = 1:n-1
        left1(i,k) = r1(i) - sum(ank1(i,1:k));
        right1(i,k+1) = r1(i) - sum(ank1(i,1:k));
    end
end
for i = 1:m
    right2(i,1) = r1(i);
    left2(i,2*n) = r1(i+1);
    for k = 1:2*n-1
        left2(i,k) = r1(i) - sum(ank2(i,1:k));
        right2(i,k+1) = left2(i,k);
    end
end
for i = 1:m
    right3(i,1) = r1(i);
    left3(i,4*n) = r1(i+1);
    for k = 1:4*n-1
        left3(i,k) = r1(i) - sum(ank3(i,1:k));
        right3(i,k+1) = left3(i,k);
    end
end


for i = 1:m
    trap(1,i) = 1/2*(ank1(i,1)*(feval(f,right1(i,1))) + ank1(i,n)*feval(f,left1(i,n)));

    for k = 2:n-1
        trap(1,i) = trap(1,i) + ank1(i,k)*(feval(f,right1(i,k))+feval(f,left1(i,k)))/2;
    end
end
for i = 1:m
    trap(2,i) = 1/2*(ank2(i,1)*(feval(f,right2(i,1))) + ank2(i,2*n)*feval(f,left2(i,2*n)));

    for k = 2:2*n-1
        trap(2,i) = trap(2,i) + ank2(i,k)*(feval(f,right2(i,k))+feval(f,left2(i,k)))/2;
    end
end
for i = 1:m
    trap(3,i) = 1/2*(ank3(i,1)*(feval(f,right3(i,1))) + ank3(i,4*n)*feval(f,left3(i,4*n)));

    for k = 2:4*n-1
        trap(3,i) = trap(3,i) + ank3(i,k)*(feval(f,right3(i,k))+feval(f,left3(i,k)))/2;
    end
end
simp1 = (4*sum(trap(2,1:m)) - sum(trap(1,1:m)))/3;
simp2 = (4*sum(trap(3,1:m)) - sum(trap(2,1:m)))/3;
boole = (16*simp2 - simp1)/15