In general, such a brute-force solution is used only if a direct one is not available. This is not the case.
For t=0 the number of the roots of the equation is either <= 4 or infinity (so 6 is out of the question).

It remains t=1 which was solved in your previous post [or, just form the 4 quadratic equations (obtained from abs(u) = +- u) and find the conditions to have integer roots]. The case t=2 is similar.

I see two problems in your worksheet.

1. The distribution corresponding to the RV MIXTURE is not correctly defined. It should be:

pdf := (p,x) -> p*PDF(U1, x)*Heaviside(x+3) + (1-p)*PDF(U2, x)*Heaviside(x-1);
cdf := (p,x) -> int(pdf(p,y), y=-infinity..x);

MIXTURE := p -> RandomVariable( Distribution( PDF = (x -> pdf(p,x)), CDF = (x -> cdf(p,x)), Conditions=[p >=0, p <=1] )):

2. It seems that Sample works only for distributions for which the support is an interval and the PDF is differentiable here.
The support of  MIXTURE is (-3,-1) union (1,3).
[I don't know how to declare such a support, or even if it's possible].
So,  Sample(MIXTURE(1/3), 10);  will fail, but

Sample(MIXTURE(1/3), 10, method=[envelope, range=1..3]);  # works.

Probably it's possible to use somehow method=custom.

I have set interface(typesetting=standard), because typesetting=extended  exhibits (here too) its weakness, not related to our problem.

1. Let's take the following version of your example:

restart;
interface(typesetting=standard):
#interface(typesetting=extended):
j0:=1; n0:=0;
V:=Vector(n0):
for j from 1 to 5 do
  try
    print('j'=j);
    V(j):=j/(j-j0)
  catch:
    print("Ante",lastexception); 
    V(j):=17; # clears lastexception apparently
    print("Post",lastexception);
  end try
end do:
V, lastexception;

It is equivalent to the original one.
V(j):=17  will clear lastexception iff j0 > n0. So, it seems that the extension of a rtable is inside of some try ... end try and lastexception is reset, but only locally.
Outside try .. end try,  lastexception has always the correct value.

2. For the second example, a workaround would be:

restart;
interface(typesetting=standard);
lastexception;
5/0;
lastexception;

5/5;
lastexception;
x:=47;
lastexception;
table();
lastexception;
[1,2,3];
lastexception;
# rtable();
try  rtable(); catch : end try ;
lastexception; # Now NOT cleared!

The arithmetic with infinity is very delicate.
The automatic simplification must be lite, so there are compromises.
When an expression is suspected to contain infinity,  some kind of simplification is recommended.
It is strange that even is must be sometimes preceded by simplify:

eq := Pi*infinity + I*infinity = infinity + I*infinity:   #  sqrt(2) too
is(eq), is(simplify(eq));

                          false, true

(Maybe this is a consequence of the fact that a constant for complex infinity is missing in Maple.)

Edit. I add some other stange results about arithmetic with oo:

limit( (x^2+1) + x*I, x=infinity );
                            infinity
limit( (x^2+1) + x^2*I, x=infinity );
                        (1 + I) infinity
limit( exp(x) + exp(x+1)*I, x=infinity );
                   -infinity I (-exp(1) + I)
simplify(%);
                   -infinity I (-exp(1) + I)
is(%=infinity+I*infinity);
                             false

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;