@emendes 

interface(rtablesize=infinity):

aux:=[-x1*x2+u0,x1*x1-2*x2,u1,u2,u3,u4,u5,u6,u7,u8,u9,u10];

f:=Vector[row](aux);

@taro  Christian's code gives a more complete factorization. But the same can be made more simple and understandable code:

e:= g^((2*(-sigma+k+1))/(-1+sigma))-tau^2:

combine(factor(expand(e)));

or

normal(combine(factor(expand(e))));

without the factor  -1  ahead.

@Carl Love   Thanks. Take a look at my second method.

@Carl Love  I think the gold circle is the set of points of tangency of black circles.

@shani2775  Use the code below

restart;

u[0] := (4/3)*c^2*cos((1/4)*x)^2; alpha := 2;

iteration := 3;

for k from 0 while k <= iteration do

u[s] := eval(u[k], t = xi); u[k+1] := simplify(u[k]-(int(diff(u[s], [`$`(xi, alpha)])+diff(u[s]*u[s], x)+diff(u[s]*u[s], x, x, x), xi = 0 .. t)))

end do;

@Carl Love   From the original within Maple we can see

f := unapply((1), x, y)

@Markiyan Hirnyk  The code of  LSSF  was written and works in 1D input, but does not work in 2D. Replace  op(Ind minus {Ind[i]})=~0  by  op(Ind minus {Ind[i]}=~0)  in the code of the procedure. Now there no any problems.

@Carl Love  You wrote "As far as I recall, a true return from is is trustworthy."  However, any Maple's bug can be rewritten in terms of  is .

Example:

is(int(1/(x^4+4), x=0..1) < 0.2);

evalf(Int(1/(x^4+4),x=0..1));

                                                              true

                                                      0.2389834593

@Earl  You wrote "But, why are your values for DthetaZero neither the under/over values I found nor the exact value which Carl Love found above?"  But this value  

1.034385*Pi;

                                                 3.249616317

is also not exact (approximate only) and the estimate is true  

               DthetaZero1=3.2496<3.249616317<3.2497=DthetaZero2

@Preben Alsholm  Thank you for the detailed analysis of the issue.

@vv  Thank you for a workaround. There seems to be incorrectly evaluated the any integral of the series  int(1/(x^4+p), x=0..1)  for  p>1 .

Examples:

seq([evalf(int(1/(x^4+p), x=0..1)),evalf(Int(1/(x^4+p), x=0..1))], p=0.5..5, 0.5);

[1.589368880, 1.589368880], [.8669729873, .8669729873], [.3059112583, .6005615062], [.2368336857, .4603941667], [.1932460853, .3735979830], [.1632003072, .3144676426], [.1412235518, .2715577471], [.1244471178, .2389834593], [.1112198369, .2134044587], [.1005233062, .1927820248]

@Carl Love  Thanks. Vote up.

@Carl Love  I wonder why the first variant for symbolic case works, and the second no:

a:=n->piecewise(n::odd, 1/n, -1/n^2):

a1:=n->`if`(n::odd, 1/n, -1/n^2):

Examples:

a(2*n-1) assuming n::integer;

a1(2*n-1) assuming n::integer;

                                                      1/(2*n-1)

                                                  -1/((2*n-1)^2)

@Earl  1. Yes it does.

2. Your calculation for the case of a solid top is right. The same result was obtained using a triple integral:

restart;

f:=unapply(5*sin(z)*exp(-2*(Pi-z)/3), z):

F:=unapply(sqrt(f(z)^2-x^2), x,z):

int((x^2+y^2), [y=-F(x,z)..F(x,z), x=-f(z)..f(z), z=0..Pi]);  # Calculation by a triple integral

int((1/2)*Pi*f(z)^2*(f(z)^2), z = 0 .. Pi);  # Your calculation

evalf(%);

@toandhsp 

map(t->`if`(nops(t)=3,[P(t[1]),P(t[2]),t[3]],P~(t)), [[-4*x-3*y-44 = 0, -4*x+3*y-14 = 0, [x = -29/4, y = -5]], [-4*x-3*y-44 = 0, -3*x-4*y-48 = 0, [x = -32/7, y = -60/7]], [-4*x-3*y-44 = 0, 3*x+4*y-2 = 0, [x = -26, y = 20]], [-4*x-3*y-44=0,4*x+3*y-6=0]]);