@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]]);

@Markiyan Hirnyk  You are right. Vote up. Your approach solves the problem for  n=7  for 1-2 minutes.

restart;

with(GraphTheory):

L := [NonIsomorphicGraphs(7, output = graphs, outputform = adjacency)]:

CodeTools[Usage](plots[display](Matrix(116,9, DrawGraph~(Graph~(L),style= circle))));

@Markiyan Hirnyk  I used your code for  n=5  and got a strange result (in Maple 2016.1):

restart;

with(GraphTheory):

L := [NonIsomorphicGraphs(5, output = graphs, outputform = adjacency)]:

seq(DrawGraph(Graph(c)), c = L);

                                    [Length of output exceeds limit of 1000000]

@Carl Love  Vote up. It is interesting that in earlier versions of Maple  IsIsomorphic  works properly for the case  n=4  also. I'm only slightly corrected your procedure, that there were no empty rows in your Matrix (for examples for cases  n=2 or n=4):

AllGraphs1:= proc(n::posint)

uses C= combinat, GT= GraphTheory, LT= ListTools, P= plots;

local 

     Gs:= map2(GT:-Graph, n, [C:-powerset(C:-choose({$1..n},2))[]]),

     NP:= plot(()-> 0, -1..1, color= white, axes= none),

     ng, ns;   

     Gs:= map(`?[]`, [LT:-Categorize(GT:-IsIsomorphic, Gs)], [1]);

     ng:= nops(Gs);

     ns:= ceil(evalf(sqrt(ng)));

     P:-display(Matrix(`if`(ng<=ns^2-ns,op([ns-1,ns]),ns), [GT:-DrawGraph~(Gs, style= circle)[], NP $ `if`(ng<=ns^2-ns,ns^2-ng-ns,ns^2-ng)]))    

end proc:

Example for n=4 (in Maple 16):

AllGraphs1(4);

 Addition.  See this link for the issue.

For some reason, this procedure in Maple 16  (for example for n=6)  is 2.5 times faster than in Maple 2016.1.

In Maple 16:

CodeTools[Usage](AllGraphs1(6));

memory used=6.11GiB, alloc change=0 bytes, cpu time=105.25s, real time=105.45s

In Maple 2016.1:

CodeTools[Usage](AllGraphs1(6));

memory used=11.29GiB, alloc change=96.00MiB, cpu time=4.65m, real time=4.77m, gc time=18.95s