a=5;
b=3;

tol=1e-4;

n=3;
innerP=P_inner(a,b,n);
outerP=P_outer(a,b,n);
aveP=(innerP + outerP)/2;
fprintf('%3d  %18.12f \n', n,aveP);

while (outerP - innerP) > 2*tol
    n=n+1;
    innerP=P_inner(a,b,n);
    outerP=P_outer(a,b,n);
    aveP=(innerP + outerP)/2;
    fprintf('%3d  %18.12f \n', n,aveP);
end