I was able to do it in under a second by using simplify with side relations (see ?simplify,siderels ). I refer to it here by its syntax skeleton simplify(..., [...]). Unlike regular simplify---which may use heuristic ideas of what is simpler---simplify(..., [...]) uses a pure polynomial algebraic algorithmic approach. 

Replacing a common subexpression with a new variable is a type of side relation. Here are the substitutions that you wanted to apply:

There are a number of factors that can make a vast difference in the amount of time and memory that the algorithm uses. These factors are referred in the help page as "monomial orderings", but the help page fails to mention just how important these are. The naive approach to your problem would be

simplify(kappa, [newpar]);

I let this run several minutes until I was nearly out of memory (it was using about 4Gig) before I killed it. Then I tried applying the side relations individually from most complicated to least (this last point is very important), like this:

newkappa:= kappa:
for k from nops([newpar]) by -1 to 1 do
     newkappa:= simplify(newkappa, [newpar[k]])
end do:
newkappa;

This works nearly instantaneously (0.140 secs) and achieves a great simplification: The only variables in the result are the s's. I also did the C=4, K=7 case in 0.656 secs. So, there you go: I took it from several days of computation to under a second.

Here's the modified worksheet:

sC.mw

In the line before the else, there appears #Iota.... The # comments out the rest of the line, so it unbalances the parentheses.

Your erroneous result seems to be related to a bad simplification of factorials. After evaluating at x=1 and r=1, your summand has the subexpression

ex:= (1+k)!/k!^2/(2-k)!:

which simplifies as

simplify(ex);

Obviously this is wrong for any integer k. A workaround is assume(k::nonnegint).

simplify(ex) assuming k::nonnegint;

To get rid of all of the RootOfs at once, you could do

evalindets(newpar, RootOf, x-> allvalues(x)[1]);

First method (my preference):  second:= eval(t, test[2]);

Second method: second:= rhs(test[2][]) or rhs(test[2][1]) or op([2,1,2], test)

Unfortunately, there is no analog to the first method for selecting the first value; it only works on expressions of the form (variable = value).

The trick is to first use command combine to bring the A into the sum, and then use algsubs(A*W[n](t)=5, ...). Finally, use expand to bring constants back out of the sums.

A*B*Sum(W[n](t)*q[n](z), n= 1..infinity);
                     /infinity               \
                     | -----                 |
                     |  \                    |
                     |   )                   |
                 A B |  /     W[n](t) q[n](z)|
                     | -----                 |
                     \ n = 1                 /
combine(%);
                  infinity                   
                   -----                     
                    \                        
                     )                       
                    /     A B W[n](t) q[n](z)
                   -----                     
                   n = 1                     
algsubs(A*W[n](t)=5, %);
                      infinity           
                       -----             
                        \                
                         )               
                        /     5 B q[n](z)
                       -----             
                       n = 1             
expand(%);
                         /infinity       \
                         | -----         |
                         |  \            |
                         |   )           |
                     5 B |  /     q[n](z)|
                         | -----         |
                         \ n = 1         /

Let me know if that's enough explanation, or if you need a more-detailed example.

The default solution method rkf45 fails miserably at this task. (I used x(4)=10, u(4)=10, z(4)=10.) Using method= dverk78 gives about 5 digits of accuracy. I am trying some other methods.

restart:
sys:= {diff(u(t),t) = 4*u(t), diff(x(t),t) = 6*x(t)+2, diff(z(t),t) = x(t)+u(t)}:
Sol1:= dsolve(sys union {u(4)=10, x(4)=10, z(4)=10}, numeric):
Sol1(0);
         [t = 0., u(t) = HFloat(1.11364764942213e-6), 
           x(t) = HFloat(-0.3333333330052356), 
           z(t) = HFloat(7.1111113895777045)]

Sol2:= dsolve(sys union eval({u(0)=u(t), x(0)=x(t), z(0)=z(t)}, Sol1(0)), numeric):
Sol2(4);
          [t = 4., u(t) = HFloat(9.742963826802734), 
            x(t) = HFloat(7.58924709847216), 
            z(t) = HFloat(9.533948806446041)]

Sol3:= dsolve(sys union {u(4)=10, x(4)=10, z(4)=10}, numeric, method= dverk78):
Sol3(0);
     [t = HFloat(0.0), u(t) = HFloat(1.12536549286537e-6), 
       x(t) = HFloat(-0.33333333294284945), 
       z(t) = HFloat(7.111111392517568)]

Sol4:= dsolve(

     sys union eval({u(0)=u(t), x(0)=x(t), z(0)=z(t)}, Sol3(0)),

     numeric, method= dverk78

):
Sol4(4);
      [t = HFloat(4.0), u(t) = HFloat(9.999698295626974), 
        x(t) = HFloat(10.000134201687256), 
        z(t) = HFloat(9.999946940854626)]

In terms of elegance and simplicity of the algorithm used, it'd be hard to beat rem:

f:= x^5+a*x^4+b*x^3+c*x^2+d*x+e:
rem(f, x^3-1, x);
                          2                    
                 (c + 1) x  + (a + d) x + b + e

@casperyc 

It really surprises me that this would happen because I simplified it.

Here's a simple example that shows why not to simplify. These are the same type of expressions as occur in your Matrices (the ones in your Matrices are longer).

simplify(1+(a+b+c+d)*(e+f+g+h)*(i+j+k+l));

 a e i + a e j + a e k + a e l + a f i + a f j + a f k + a f l
 + a g i + a g j + a g k + a g l + a h i + a h j + a h k
 + a h l + b e i + b e j + b e k + b e l + b f i + b f j
 + b f k + b f l + b g i + b g j + b g k + b g l + b h i
 + b h j + b h k + b h l + c e i + c e j + c e k + c e l
 + c f i + c f j + c f k + c f l + c g i + c g j + c g k
 + c g l + c h i + c h j + c h k + c h l + d e i + d e j
 + d e k + d e l + d f i + d f j + d f k + d f l + d g i
 + d g j + d g k + d g l + d h i + d h j + d h k + d h l + 1

Before simplification, there are 2 multiplications and 10 additions; after, there are 128 multiplications and 64 additions. So the complexity increased by a factor of 16. One way to guard against this is to include the size modifier to simplify.

simplify(1+(a+b+c+d)*(e+f+g+h)*(i+j+k+l), size);
      1 + (a + b + c + d) (e + f + g + h) (i + j + k + l)

But, based on an extremely brief viewing of some of the entries in your Matrices, I guess that simplify will provide no benefit in your case, even with the size option.

The functional way:

f:= x-> exp(-exp(x^2)):
plot(z-> Int(f, 0..z), -3..3);

Note the capital I in Int and the lack of a variable appearing before the ranges.

I assume that you intended for g to be considered a function of x, y, z. In that case, you should change Diff(g, x) to D(g)(x) or Diff(g(x,y,z), x), and likewise for y and z.

When you simplify the expressions in these matrices, they become incredibly long---essentially, they are expanded as polynomials in the ~ 10-20 variables. The only change I made was to remove the line kappa:= simplify(kappa). With this change, I was able to eval the K=8 case (the 254x20) in 0.25 seconds, and I did the K=10 (1022 x 24) case in 2.5 seconds and the K=11 (2046 x 26) in 8.5 seconds. The mykappa itself also runs much, much faster without the simplify. So, there, I gave you a factor of about 1000 improvement.

Not simplifying should also give you greater accuracy in the floating-point evaluation---fewer operations means less round-off error.

With some research in Wikipedia, I found an algorithm for counting the Eulerian circuits (not paths) in directed graphs (digraphs). It's called the BEST algorithm, BEST being the initials of the four authors' surnames. A digraph has an Eulerian circuit iff for every vertex v, indegree(v) = outdegree(v). When these are equal, we will just call it deg(v). Denote by N(G) the number of spanning trees of G. Then the number of Eulerian circuits is

N(G)*mul((deg(v)-1)!, v in G)

(For the graph in question, clearly all the factorials are 1.) Maple has a command for counting the spanning trees: GraphTheory:-NumberOfSpanningTrees. Applying the BEST algorithm in this case, we get 

restart:
macro(GT= GraphTheory):
G1:= GT:-Graph( {[C,G], [G,C], [T,G], [G,T], [A,T], [C,A], [T,C]}):
GT:-NumberOfSpanningTrees(G1)*mul((v/2-1)!, v in GT:-DegreeSequence(G1));
                               3

So, there are three Eulerian circuits.

(The division by 2 in the factorial is because DegreeSequence adds together the InDegree and OutDegree.)

See this answer from last week: http://www.mapleprimes.com/questions/150074-Invalid-Left-Hand-Side-In-Assignment-Problem#comment150076

There is no bug. Option inplace, the default, means that LinearSolve does not return a value directly; instead, it alters the input matrix.

(Dammit, the worksheet uploader is broken right now! I have to upload this worksheet in plaintext without the output. You'll have to download the worksheet to see the output.)

restart:
macro(LA= LinearAlgebra, LAM= LinearAlgebra:-Modular):
# Generate a random example linear system AX=B. The purpose of these while loops is to make sure that
# the example is not trivial.
det:= 0:
while det = 0 do
     A:= LAM:-Create(2, 3, 3, random, integer[]);
     det:= LAM:-Determinant(2, A)
end do:
N:= 0:
while N=0 do
     B:= LAM:-Create(2, 3, 1, random, integer[]);
     N:= LA:-Norm(B)
end do:

Aug:= < A | B >;
              
# By default, Modular:-LinearSolve works inplace, which means that it does not return a value directly.
# Instead, it alters the input matrix.
LAM:-LinearSolve(2, Aug, 1);
Aug;
             
X:= Aug[..,-1];
          
# Verify solution.
LAM:-Multiply(2, A, X);
          
B;
              
# Contrast with solution obtained with inplace= false.
X:= LAM:-LinearSolve(2, < A | B >, 1, inplace= false);

#Download the worksheet:  LA_Modular.mw