The issue that you're experiencing is specific to RootOf's use of _Z; it does not apply to _Z in general, and it has nothing to do with the underscore. Once you understand the issue, it's very easy to avoid.

First, let me explain bound variables and free variables. These are fundamental concepts of mathematical logic; they aren't specific to Maple or even to computer languages. For example, consider the commands 

int(a*x^2, x), int(a*x^2, x= 1..2),

in the absence of any surrounding context that may change what a and x mean. In both commands, a is a free variable and x is a bound variable. This is true regardless of whether they are local or global. The key thing is that x cannot be assigned any value (other than another name) without turning the command into nonsense.

Many commands require a bound variable or variables: int, diff, limit, sum, seq, solve (usually), plot (usually), RootOf, and many others. (In cases where the main argument is allowed to be a procedure, we can think of the procedure's parameters as the bound variables.) RootOf is unique among these (unfortunately) in the way that it handles its bound variable. If you do not specify a bound variable as the 2nd argument, it will look for global _Z in the first argument. If it's not there, it's an error. If you do specify a bound variable, RootOf will always change it to global _Z. Since _Z is protected, it should never be assigned a value, and you don't need to worry about that. So, in order for your example with local _Z to work, you need to put that _Z as the second argument of RootOf. (But there's really no reason to do this because you'd might as well just use global _Z, and, regardless, it's going to change whatever you use to global _Z anyway.)

The statement given at ?_ about those variables being reserved is ridiculously overstated. For one thing, it's not possible for a local variable beginning with _ to cause any problem.

You should never declare a global variable as global unless you intend to assign a value to it within the procedure containing the declaration.

You need to replace e^(-B.eta) with exp(-B*eta).

The number of boundary conditions required is the differential order, 6, plus the number of unspecified constants, e.

Maple does not "struggle with integer factorizations"; you just happened to stumble upon a particular number that exposes a bug. In the worksheet below, I randomly generate 1000 composites, each a product of a 9-digit prime and an 11-digit prime (like your example), and factor them all in 7 seconds.

Download ifactor.mw

When mod is used on a polynomial, such as t, it's mapped over the coefficients; so, t mod 10 becomes just t in the 2nd plot. This happens before plot even sees its arguments.

In the first plot, the special evaluation rule of seq prevents the simplification of t mod 10 to t.

This will do that job for any list:

ListTools:-Collect(L)

You are essentially creating an NxN matrix a. Since you don't explicitly make it a matrix, it becomes a table. A table is not so bad, but a matrix would be better. The source of your error is that the diagonal entries a[i,i] must be computed after the other elements. In the worksheet below, I show two methods for this. In deriving the second method, I completely eliminated the need for the polynomial and its derivative.

Download BurgersMatrix.mw

Your p can be factored as p = (6000 - w^4)*(z*g__u + f*g__z)/144000, so p=0 implies that w doesn't depend on any of the other variables. w = 6000^(1/4), which simplifies to 2*375^(1/4).

y is a name, not a string; "y" is a string, not a name. They are not the same thing.

fx:= 3*x^2 - 9*y:
str:= String(fx);
                       str := "3*x^2-9*y"
has(fx, y);
                              true

evalb(SearchText("y", str) <> 0);
                              true

Why are you converting mathematical expressions to strings? I suspect that there is an easier way to achieve your ultimate goal.

There's no technical reason why those errors are untrappable; the reason is purely philosophical, and it's stated on help page ?try. Here is an excerpt from the 10th paragraph of Description:

• The following exceptions are untrappable:
  - Computation timed out (this can only be caught by timelimit, which raises a "time expired" exception, which can be caught).
  - Computation interrupted (user pressed Ctrl+C, Break, or equivalent).
  - Internal system error (which indicates a bug in Maple itself).
  - ASSERT or return type or local variable type assertion failure (assertion failures are not trappable because they indicate a coding error, not an algorithmic failure).
  - Stack overflow (when that happens, generally stack space is insufficient for doing anything like running cleanup code). This includes the "too many levels of recursion" exception at the Maple language level.

Pick a value of x larger than the stated singularity, say xs. Reformulate the initial conditions to be at xs rather than 0. Call dsolve again. Plot from xs to 10 (or wherever). You may encounter new singularities and need to repeat this process.

orthopoly is an old package, created before modules existed. It is stored in a table. The A:-B syntax works for table-based packages, but apparently proc() uses A; B() end proc() does not.

You have 

for {w = 30..150}
    s:= ...;

How did you come up with that syntax? Change it to 

for w from 30 to 150 do
    s[w]:= ...
end do;  # "end do" can be abbreviated "od"

I used s[w] instead of s so that the results will be saved for each w.

Error messages that are accompanied by a red dashed box result from very basic syntax errors in the command that is currently being entered. They can't possibly have anything to do with code already entered. Thus, this error has nothing to do with piecewise.

You could read the help page ?for. Also, nearly any freely available AI chat assistant will give a reasonable answer to "How do I write a for loop in Maple?"

Here is how I think chaotic IVP solutions should be plotted. Note the use of the options parameters, range, and refine. The differences between these plots and the plots shown in the paper are too categorical to possibly be explained by any accumulation of rounding errors or chaotic process. For example, the second plot below shows only positive values for x(t) whereas the corresponding plot in the paper shows only negative values for x(t).

I've tried several different numeric methods (rkf45, ck45, rosenbrock, dverk78) and different levels of error control (relerr= 5e-6, abserr= 5e-6, relerr= 5e-9), but I got the same plots.

Download DuffingEq.mw

 As always, plots look much crisper in the actual worksheet than they do when the worksheet is uploaded to MaplePrimes.

The command you need is convert(f(x), fullparfrac, x), possibly with added options. See the help page ?convert,fullparfrac.

Maple has a "sum over all roots of a polynomial" feature that usually makes seeing the individual factors unnecessary. However, if that's not useful for you, try including the factor option.

You get much better simplifications by calling dsolve for each numeric instantiation rather than getting a generic solution and then substituting parameter values,

Download DataFRame.mw