The function f is strictly convex (because x -> x^2 is strictly convex and the other 4 terms are convex).
So f has at most 2 real solutions (and you want 6 integer ones).

msolve(x^3+y = 123, 124);
        {x = x, y = 123*x^3 + 123}

The equation is trivial because x is obviously arbitrary.
Take a nontrivial one:

msolve(x^3+y^2 = 123, 124);
   {x = 7, y = 20}, {x = 7, y = 42}, {x = 7, y = 82},  {x = 7, y = 104}, {x = 11, y = 30}, {x = 11, y = 32},
    {x = 11, y = 92}, {x = 11, y = 94}, {x = 27, y = 30}, {x = 27, y = 32}, {x = 27, y = 92}, {x = 27, y = 94},
    {x = 35, y = 20}, {x = 35, y = 42}, {x = 35, y = 82}, {x = 35, y = 104}, {x = 39, y = 18}, {x = 39, y = 44},
    {x = 39, y = 80}, {x = 39, y = 106}, {x = 43, y = 10}, {x = 43, y = 52}, {x = 43, y = 72}, {x = 43, y = 114},
    {x = 51, y = 20}, {x = 51, y = 42}, {x = 51, y = 82}, {x = 51, y = 104}, {x = 55, y = 30}, {x = 55, y = 32},
    {x = 55, y = 92}, {x = 55, y = 94}, {x = 71, y = 18}, {x = 71, y = 44}, {x = 71, y = 80}, {x = 71, y = 106},
    {x = 83, y = 10}, {x = 83, y = 52}, {x = 83, y = 72}, {x = 83, y = 114}, {x = 91, y = 10}, {x = 91, y = 52},
    {x = 91, y = 72}, {x = 91, y = 114}, {x = 99, y = 0}, {x = 99, y = 62}, {x = 107, y = 18}, {x = 107, y = 44},
    {x = 107, y = 80}, {x = 107, y = 106}, {x = 119, y = 0}, {x = 119, y = 62}, {x = 123, y = 0}, {x = 123, y = 62}

Inside a loop, the two statement separators (semicolon and colon) are equivalent. Ony the terminator after end counts (semicolon in your case). See:
?for
?;
?printlevel

I don't know what you mean by "not well posed" in this context. I think that it is more important to state clearly the class of the function.

So, if u is supposed to be differentiable in (0, oo) x (0, oo) and continuous in [0, oo) x [0, oo) then the problem has a unique solution:

u(x, t) =  piecewise(x>t, exp(-(x-t)^2), 1);

Maple does not find it of course, but does half of the job, by solving the pde alone.
(bc is interpreted as D[1](u)(0+, t) = 0 for t>0).

In Maple 2019 does
(1/2)*(int(cosh(cos(2*t)), t = 0 .. 2*Pi) + int(cosh(-cos(2*t)), t = 0 .. 2*Pi));
produce 2*Pi*BesselI(0, 1)?  Normally, it should simplify cosh(-cos(2*t)) to cosh(cos(2*t)).

In Maple 2018 I have to manipulate like this:

plot3d({8 - x^2 - y^2,x^2 + y^2}, x=-2..2, y=-sqrt(4-x^2)..sqrt(4-x^2));

The plot shoul look like

obtained with

plots:-inequal([x^2+y^2 >= x, x^2+y^2 <= 1/9], x=-1/3..1/3, y=-1/3..1/3,
optionsfeasible=[color=red], optionsclosed=[color=white], scaling = constrained);

The option coords=...  fails for some plot commands.

It's a version of the classical:

MDS:=proc(A::Matrix, r::integer) 
local n,B,i,j,ip,jp,k,P,W, rng;
rng:=t -> (t-1)*r+1 .. t*r; 
n:=upperbound(A,1)/r;
for i to n do   for j to n do 
  B[i,j]:=LinearAlgebra:-SubMatrix(A, rng(i), rng(j), datatype=integer[8])
od od;
for k to n do
   W := Matrix(r*k, datatype=integer[8]); 
   P:=combinat:-choose(n,k);
   for ip in P do for jp in P do
      for i to k do for j to k do
         W[rng(i),rng(j)] := B[ip[i],jp[j]];
      od od;
      if LinearAlgebra:-Modular:-Determinant(2,W)=0 then return 'NoMDS' fi;
   od od;
od;
'IsMDS'
end proc:

A:=Matrix(   # check
 < 1,0,1,0;
   0,1,0,1;
   1,0,0,1;
   1,1,1,0>):
MDS(A,2);
A[1,1]:=0:
MDS(A,2);
                             IsMDS
                             NoMDS

If you are working in Z2 (i.e. mod 2) then a square matrix A is MDS iff A = [1]  (so the dimension is 1).
MDS matrices are useful only for fields other than GF(2).
So, the MDS procedure reduces to

MDS:=proc(A::Matrix, r::integer)
  if r<> 1 then return 'NONE' fi;
  eval(A, 0=NO)
end proc;

e_n_1b := n_1 = (-w^(2*sigma)*tau + w^sigma)*s*nsp_1/(w^sigma*tau^2 - w^(2*sigma)*tau + w^sigma - tau) + tau*(w^(sigma - 1)*tau - w^(2*sigma - 1))*s*nsp_2/(w^sigma*tau^2 - w^(2*sigma)*tau + w^sigma - tau):
lhs(e_n_1b)=map(factor@expand, rhs(e_n_1b))

You want to simplify a polynomial f in n (=7) indeterminates.
Maple has simplify(f, size)  for this (size is by default).
Of course the size is not necessarily minimal, and AFAIK such algorithm does not exist.
But if you know/hope/suspect a special convenient form for f, tell us about it.

restart;
f:=
arccos((-sin(theta3)*cos(phi3)*cos(theta2)+cos(theta3)*sin(theta2))/((
(1-cos(theta3)^2)*cos(theta2)^2+cos(theta3)^2-1)*cos(phi3)^2-2*sin(theta3)*cos(
theta2)*cos(theta3)*sin(theta2)*cos(phi3)-cos(theta2)^2*cos(theta3)^2+1)^(1/2))
*sin(theta3):

evalf[15](Int(Re(f), [theta2=0..Pi,theta3=0..Pi,phi3=0..2*Pi], method=_CubaCuhre, epsilon=1e-6));
                        62.0124733278130
evalf[15](Int(Im(f), [theta2=0..Pi,theta3=0..Pi,phi3=0..2*Pi], method=_CubaCuhre, epsilon=1e-6));
                       Float(undefined)