@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.

Your function f : R^3 --> R should be locally integrable, but it's not!
So, the Hardy−Littlewood maximal function makes no sense.

@Thomas Richard Sometimes the situation is reversed!

@sursumCorda So, you want an explicit solution instead of a parametric one (t3 being the parameter).
This could be very difficult (or even impossible) to obtain in general.

It would be nice to tell us how the equation was obtained.

Here is a (simpler) parametric one (Klimek G.&M. - Discovering curves and surfaces with Maple, Springer 1997.)
 

plot3d([(4+3.8*cos(y))*cos(x),
        (4+3.8*cos(y))*sin(x), 
        (cos(y)+sin(y)-1)*(1+sin(y))*log(1-Pi/10*y)+7.5*sin(y)],
      x=-Pi .. Pi, y=-Pi .. Pi, color=[sin(x+1)+cos(y), sin(y), 0.3],
      light=[85, 50,  0.6, 0.8, 0.2], orientation=[47, 55],
      style=patchnogrid, axes=none);