The function f defined by fun is 0-homogeneous ( f(tx,ty,tz) = f(x,y,z), for t>0 )  and symmetric, so, in order to compute the sup we may consider just f(x,y,1)  for (x, y) in [0,1]^2 \ A.

Note that 

limit(f(a^2, a,1) , a=0, right) = 3,  limit(f(a, a,1) , a=0, right) = 4*sqrt(2) - 3 = 2.65... 

So, the sup is probably indeed 3, but I wonder what numerical test could confirm it at least partially! 

@mmcdara The fact that T(S) has a power of 2 elements is a classical result (for arbitrary cardinals); as a consequence, T(S) cannot be countable. I used to prove this to my students.

The number of sets in a (sigma-)algebra should be a power of 2 i.e. 2,4,8, ...

Yes, but PDEtools:-Solve shows a nonsense instead.

@Carl Love A very professional code. Vote up!

1. You forgot to state the (optimal control?) problem.

2. ss seems to contain an algorithm for the unstated problem. Are you sure that it is correct?

3. Do you know the exact solution of your problem. Probably it can be obtained because the ODE has a simple symbolic solution.

4. NLPSolve does not work for me for nvar=3. Does it work for you as it is?

5. Have you considered the fact that NLPSolve gives in general only local minima?
 

@tomleslie See the matrix form of the command.

@Carl Love Yes, but then the results of LC and LeafCount will not be exactly the same.

@sursumCorda Yes, I see it now.

It seems that in seq(expr, 1..n),  expr is evaluated n times   but in seq(expr, n) ,  expr is evaluted n+1 times (being evaluated once befor creating the sequence).
But note that this is not actually a bug,  because the second form is not documented!

restart;
u:=10;
seq((u:=u+1), 3)
#                            u := 10
#                           12, 13, 14
u:=10;
seq((u:=u+1), 1..3)
#                            u := 10
#                           11, 12, 13

@sursumCorda 

evalf(Int(convert(x^(x^x), exp), x=1..6));

     1.102664999 * 10^36300
 

@Carl Love I obtain an error due to complex values for sample. But this works:

plot([LambertW(x), LambertW(-1,x)], x=-1..4, view=[-1..4, -3.5..1.5], colour=[red,blue], scaling=constrained, labels=[x,LambertW(x)], adaptive=15);

@Axel Vogt  mma gives:

The problem here is that LerchPhi is evaluated numerically very slowly for large arguments.
Yes, this could be seen as a bug, but there are easy workarounds. For example,

L:=LerchPhi(1/10, 1, 1/2-100000000*I);
evalf(L);   # must be interrupted

but 

evalf(convert(L,Sum)) ;  # works with any precision

@Zeineb Ok, I understand now, but you merge Maple and math notations: in piecewise I interpreted y as a second component of a point (x,y,z) in R^3, instead of a point.