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

This works: I explicitly construct the relevant polygon from the data provided by arc:

M:= op(1, plottools:-arc([0,0], 1, Pi/6..3*Pi/4)):
plots:-polygonplot(< < M[1,1] | 0 >, M, < M[-1,1] | 0 > >);

Let me know how that goes for you.

Regarding the help: Some of the options listed on the page ?plot,options apply only to the plot command itself. I think that those options are adaptive, discont, filled, resolution, sample, and smartview. And option filledregions applies only to the three plots commands contourplot, implicitplot, and listcontplot.

Assuming that the B's are nonnegative (if the B's can be negative, you'll need to explain the significance of that):

A:= [3,5,7,2,11,2,1,10,6,15]:
B:= [0.23, 0.11, 0.05, 0.05, 0.04, 0.06, 0.10, 0.20, 0.06, 0.10]:
Statistics:-Histogram(zip(`$`, A, round~(1000./`+`(B[]) *~ B)), discrete);

How about this?

restart:
with(GraphTheory):
G1:= Graph(
     [A, AT, T, CA, TC, GT, TG, C, CG, G, GC],
     {
          [C,CG], [CG,G], [G,GC],
          [GC,C], [G,GT], [GT,T],
          [T,TG], [TG,G], [A,AT], [AT,T],
          [C,CA], [CA,A], [T,TC], [TC,C]
     }
):
SetVertexPositions(
     G1,
     [
          [0,2], [1.5,2], [3,2],
          [0,1], [1.5,1], [3,1], [3.5,1],
          [0,0], [1.5,0], [3,0],
                 [1.5,-.5]
     ]
);
HighlightVertex(G1, [A,T,C,G], red);
HighlightVertex(G1, [AT,CA,TC,GT,TG,CG,GC], grey);

DrawGraph(G1);

I am guessing that by "labelled" you mean that you want the contours themselves labelled. That would perhaps be possible with some carefully placed textplots. But perhaps the following is good enough: I've placed the labels on the z-axis and some eye-tracking lines from the labels to the contours.

Contours:= [10*k $ k =-3..2]:
P1:= plots:-contourplot3d(
     x^3-y^2, x= -3..3, y= -3..2,
     axes= box, shading= zhue, filled,
     contours= Contours,
     tickmarks= [default$2, Contours],
     transparency= 0.5, thickness= 3
):
P2:= seq(
     plots:-pointplot3d(
          [[-3,2,c], [3,2,c]],
          style= line, color= gray
     ),
     c= Contours
):
plots:-display([P2,P1]);

The following command will create the graph. Having Maple draw the graph exactly the same way as you did is another matter. But Maple's graph is isomorphic to yours.
restart:
with(GraphTheory):
G1:= Graph( {[C,G], [G,C], [T,G], [G,T], [A,T], [C,A], [T,C]});
   
DrawGraph(G1);

graph.mw