The following procedure computes the sequence. However, I don't get the values that you propose. For example, how do you get 2, not 4, as the second entry?

A:= proc(n::posint)
option remember;
local A0; #A(n-1)
   if n=1 then return 1 fi;
   A0:= thisproc(n-1);
   `if`(A0::even, `if`((A0+1)::prime, A0+1, n+2), A0+3)
end proc
:
seq(A(n), n= 1..20);
1, 4, 5, 8, 10, 11, 14, 16, 17, 20, 22, 23, 26, 28, 29, 32, 34, 36, 37, 40

The name Zeta is already defined in Maple (for the Riemann-Hurwitz zeta function), so it's dangerous to use it as a variable. You should replace it with zeta, which'll prettyprint identically to Zeta as the lowercase.

Once you make that change, then all you need to do is

evalc(%)

The evalc automatically assumes that all variables are real, thus avoiding some of the complexity of the other Answers. But it won't assume that Zeta is real because that's already defined and is thus not a free variable.
 

If you'd prefer to continue using Zeta, you could do

local Zeta;

and it'll now be a totally free variable. If perchance you'd still like to use the Riemann-Hurwitz zeta function, you can call it with :-Zeta(...), and now the names Zeta and :-Zeta will not have any link to one another.

You can use lhs and rhs to operate on specific sides of equations. You can use rcurry to pass trailing arguments to operators. So, what you ask for can been done as

(lhs = rcurry(int, T[1])@rhs)(%)

I can't say whether it's "the right thing to do."

It wasn't removed. It was branched into a separate Question (entitled "Algorithms") for which your Comment became an Answer. The branching wasn't done by me, although I do agree with the editorial decision to do so. I converted the branch from a Post to a Question, converted your Comment to an Answer, and voted it up. Any Moderator can do these things; it wasn't VV.

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.