function y = closed1dK(f,a,b,c,d,n)

% f is the 2d monotone function in question
% n is the number of partitions
% m is the degree of precision for the Romberg, it is set in line 20
% a and b are our limits, wrt x, c and d wrt y
% Rom1 will return a finite result if the singularity is near a, Rom2 will
% return a result if the singularity is near b.
warning off all;

m = 3;
for k = 1:n
    
    ankx(1,k) = (b-a)*2*k/(n^2 + n);
    anky(1,k) = 2*(d-c)*k/(n^2 + n);
    %anky(1,k) = 6*(b-a)*k^2/(n*(n+1)*(2*n+1));
    %anky(1,k) = 4*(b-a)*k^3/(n^2*(n+1)^2);
    xright1(1,k) = a + sum(ankx(1,1:k));
    xright1(2,2*k) = a + sum(ankx(1,1:k));
    xright1(3,4*k) = a + sum(ankx(1,1:k));
    
    yright1(1,k) = c + sum(anky(1,1:k));
    yright1(2,2*k) = c + sum(anky(1,1:k));
    yright1(3,4*k) = c + sum(anky(1,1:k));

    

end
xright1(2,1) = a + ankx(1,1)/2;
xright1(3,1) = a + ankx(1,1)/4;

yright1(2,1) = c + anky(1,1)/2;
yright1(3,1) = c + anky(1,1)/4;

for k = 1:n-1
    xright1(2,2*k+1) = (xright1(1,k)+xright1(1,k+1))/2;
    yright1(2,2*k+1) = (yright1(1,k)+yright1(1,k+1))/2;
end

for k = 1:2*n-1
    xright1(3,2*k) = xright1(2,k);
    xright1(3,2*k+1) = 1/2*(xright1(2,k)+xright1(2,k+1));
 
    yright1(3,2*k) = yright1(2,k);
    yright1(3,2*k+1) = 1/2*(yright1(2,k)+yright1(2,k+1));

end


    

ankx(2,1) = ankx(1,1)/2;
ankx(3,1) = ankx(1,1)/4;

anky(2,1) = anky(1,1)/2;
anky(3,1) = anky(1,1)/4;

for k = 2:2*n
    ankx(2,k) = xright1(2,k) - xright1(2,k-1);
    anky(2,k) = yright1(2,k) - yright1(2,k-1);
end

for k = 2:4*n
    ankx(3,k) = xright1(3,k) - xright1(3,k-1);
    anky(3,k) = yright1(3,k) - yright1(3,k-1);
end
%{
for i = 1:3
   
    for k = 2:4*n
    
            xleft1(i,k) = xright1(i,k-1);
            yleft1(i,k) = yright1(i,k-1);
    end
    
end
%}
Rom1 = 0; Rom2 = 0; Rom3 = 0;

Rom1 = ankx(1,1)*anky(1,1)*feval(f,xright1(1,1),yright1(1,1))/4;


if m <= 2
        for k = 2:n
            for i = 2:n
                Rom1 = Rom1 + ankx(1,k)*anky(1,i)*(feval(f,xright1(1,k),yright1(1,i)) + feval(f,xright1(1,k-1),yright1(1,i-1)) + feval(f,xright1(1,k-1),yright1(1,i)) + feval(f,xright1(1,k-1),yright1(1,i)))/4;
            end
            
        
        end

        Rom1
    return
   
end



S1(1,1) = ankx(1,1)*anky(1,1)*feval(f,xright1(1,1),yright1(1,1))/4; 
S1(1,2) = ankx(2,1)*anky(2,1)*feval(f,xright1(2,1),yright1(2,1))/4; 
S1(1,3) = ankx(3,1)*anky(3,1)*feval(f,xright1(3,1),yright1(3,1))/2; 

for i = 1:m
    for k = 2:n
        S1(1,i) = S1(1,i) + ankx(i,1)*anky(i,k)*(feval(f,xright1(i,1),yright1(i,k)) + feval(f,xright1(i,1),yright1(i,k-1)))/4;
        S1(1,i) = S1(1,i) + ankx(i,k)*anky(i,1)*(feval(f,xright1(i,k),yright1(i,1)) + feval(f,xright1(i,k-1),yright1(i,1)))/4;
    end
end
%%%%%%%%%%%%%%%%%%
if m >= 2
    for i = 1:m
        for v = 2:2^(i-1)*n
        
            for k = 2:2^(i-1)*n
            S1(1,i) = S1(1,i) + ankx(i,k)*anky(i,v)*(feval(f,xright1(i,k),yright1(i,v)) + feval(f,xright1(i,k-1),yright1(i,v)) + feval(f,xright1(i,k),yright1(i,v-1)) + feval(f,xright1(i,k-1),yright1(i,v-1)))/4;
        
           
            end
        end
    end
        
    
        Rom1 = (4*S1(1,2) - S1(1,1))/3;
        T1 = Rom1;
        T2 = (4*S1(1,3) - S1(1,2))/3;
        Rom2 = (16*T2 - T1)/15;

end
Rom2