% itzung S. 174
g=inline(...
    '1./(1-2*sqrt((1/4+ep)*(1/4-ep))*cos(phi)-cos(theta)/2)/pi^2',...
    'phi','theta','ep');
f=inline('dblquad(g,0,pi,0,pi,1e-13,@quadl,ep)-2','ep','g');
epsilon=fzero(f,[0.06,0.07],optimset('TolX',1e-16),g)