Instead of 

pdsolve(PDE, IBC, ...)

use 

pdsolve(PDE, {IBC}, ...)

This will correct that error that you're currently facing. I'm not necessarily saying that this will get you to the final result.

There is no need to fully factorize the numbers. If the part that remains after factoring out all 2s, 3s, 5s, 7s, and 11s is greater than 1, then the number must have a prime factor greater than 11. Thus, the following code is much faster than the other Answers that call ifactor or ifactors:

QuoP:= proc(n,p)
local q:= n, q1;
   do if irem(q,p,'q1')=0 then q:= q1 else return q fi od
end proc
:  
QuoPs:= proc(n)
local q:= n, p;
   for p in [2, 3, 5, 7, 11] while q <> 1 do q:= QuoP(q,p) od;
   q
end proc
:
Check:= proc(n)
local k;
   for k from n to n+3 do if QuoPs(k) <> 1 then return true fi od;
   false
end proc
:

Using certain reasonably simple yet also reasonably general assumptions, that expression simplifies to 0. By eye, I see that a >= 0, b >= 0, z >= 0 are sufficient conditions. These are not necessarily necessary conditions, just an example of what I see immediately as sufficient. I think that the point being made by pdetest (a point supported by its documentation) is that the pde solution is only guaranteed to be valid under certain assumptions.

The local command can be used at the top level. So, all you need to do is

local gamma;

Your transition matrix P is a patterned tridiagonal or band matrix. It can be constructed with a single-line Matrix command using the scan= band option.

I can see no practical value to storing the positive integer powers of this matrix, so I haven't done that. You might as well simply recompute the powers as needed.
 

Download TransitionMatrix.mw

I think that you have enough sophistication and experience with a variety of languages^[*1] to realize that any statement of the form "ASCII string X and ASCII string Y mean exactly the same thing in computer language L, which is abstract mathematical object M" is ontologically problematic. The level of evaluation (parser, automatic simplifier, etc.) needs to be considered. For your example, lprint (among numerous possible evaluations) reveals that they're not the same thing at the default level of evaluation. 

The only form of derivative for initial and boundary conditions containing derivatives that's consistently supported by the documentation is D. If some component accepts another form, that's merely grace (to use Preben's word) or a bonus (to use Vv's word). So it's only at the highest level of evaluation---inside the mind of the programmer---that they could be considered equal.

^[*1]Judging by your 752 postings on MaplePrimes, every word of which I've read.

If LL is your list of lists, then the list of the averages of the sublists can computed by (add/nops)~(LL). This works for both symbolic and numeric data. 

You can't meaningfully assign values to both a name and an indexed form of the same name. In this case, you've done this with r and r[1] and k and k[1]. You can avoid this problem while still having the prettyprinted form of the index appear as a subscript by using r__1. If you don't care about whether the index appears as a subscript, you can use r and r1.

There is no problem with assigning to multiple indexed forms of the same main name as long as you don't assign to the main stem. For example, you can assign to r[0] and r[1] as long as you don't assign to r.

The command that you're trying to use, ExpandSteps, is only for polynomial expansion, not for matrix arithmetic. The command that applies to your situation is Student:-LinearAlgebra:-GaussJordanEliminationTutor; however, its use is limited to matrices at most 5 x 5.

I see that your matrix has 9-digit exact integers and exact imaginary values. So what's the point of seeing the steps? Do you suspect an error in Maple's computation?

Maple handles regular expressions. See these help pages:

• ?StringTools,Regular_Expressions
• ?StringTools,RegMatch
• ?StringTools,RegSplit
• ?StringTools,RegSub
• ?StringTools,RegSubs

All that VV said is true. I'd just like to add that if you want symbolic variables to be treated like would-be integers until they have values assigned, then use irem instead of mod.

@tomleslie Thanks for testing. The results are not surprising given that I used an array half the size of yours but a more complicated indexing scheme, and otherwise our codes are remarkably similar. Here's an update---much terser but not significantly faster---that probably satisfies the OP's titular no-loop request (thus I made this an Answer):

EratosthenesSieve:= proc(N::posint)
description "Returns list of primes <= N";
local p, k, P:= rtable(1..iquo(N-1,2), k-> 2*k+1);
   if N < 5 then return remove(`>`, [2, 3], N) fi;
   P[[seq(seq(iquo(k-1, 2), k= p^2..N, 2*p), p= thisproc(isqrt(N))[2..])]]:= ();
   [2, seq(P)]
end proc:

It is useful---and often essential---to avoid repetition when generating combinatorial objects. It can, for example, reduce the time complexity from O(n^2) to O(n), which can be the difference between doability and nondoability.

Here is a generalization of your problem: Let B be the base set from which the integers are drawn, e the exponent, and n the number of integers to draw (with replacement, but regardless of order). Then your problem can be done by:

B:= [$0..5]:  e:= 2:  n:= 2:
[seq(add(B[[seq(C)]]^~e), C= Iterator:-MultiCombination([n$nops(B)], n))];

Note that I used no sets, only lists, yet there are no repetitions. Thus, each desired object was generated only once. The Iterator package is generally very good about that.

Just change the last line to

Maximize(simplify(S), {a^2+b^2+c^2 = 1}) assuming real;

This is needed to get rid of the absolute-value / complex-modulus operators, which are messing up the differentiability of the objective function.

I can confirm using your worksheet that I get the same syntax error as you do. I don't know why. But you can do this:

q1:= (n1, n2)-> `if`(n1 < 0 or n2 < 0 or abs(n1-n2) <> 1, 0, sqrt(max(n1,n2)/2));

That gives me no syntax error.