Compact:

seq(((V^%T . W)[i,i] <= (X^%T . Y)[i,i], i=1..S);

Of course, this is not very efficient because it calculates the same matrix product over and over again.

Somewhat better:

A:= V^%T . W; B:= X^%T . Y;
seq(A[i,i] <= B[i,i], i=1..S);

Better yet: only calculate the products you need.

with(LinearAlgebra):
seq(Column(V,i) . Column(W,i) <= Column(X,i) . Column(Y,i), i = 1 .. S);

You could try this:

> subs(a2=A2,A1 . a2);

                       [  E[0] d[ba]                ]
                       [- ----------       delta    ]
                       [    ℏ                  ]
                 1/2 I [                            ]
                       [                  E[0] d[ab]]
                       [   -delta       - ----------]
                       [                    ℏ  ]

The point being that the result of subs is not evaluated before being displayed.  
However, when evaluated again, the 1/2*I will be multiplied in to the Matrix.

pdsolve({(-x-y)*diff(u(x,y),x) - y*diff(u(x,y),y)+u(x,y)=0,u(0,y)=g(y)},u(x,y),method=characteristic);

On the curve x = x(s), y=y(s), you will have du/ds = (dx/ds) u_x + (dy/ds) u_y.  Therefore if the PDE is
a(x,y) u_x + b(x,y) u_y + c u = 0, and the characteristic curve is a solution of dx/ds = a(x,y), dy/ds = b(x,y),
along that characteristic curve you have du/ds + c u = 0.  Now if x(0) = 0 and your initial condition is u(0, y) = g(y),
this gives you u(x(s), y(s)) = u(x(0), y(0)) exp(-c s) = g(y(0)) exp(-c s).

You might try DrawGraph in the GraphTheory package, with the option style=tree

What exactly do you want the C code to do?  solve may return one or more symbolic expressions for the solution.  I don't think you expect your C program to do this.  Do you have a numerical solving procedure available in some C library that you would like this program to call?  Do you want to write your own procedure C to solve the system numerically?  Or do you want your C program to evaluate the symbolic expressions that Maple's solve returns for your system (when given symbolic parameters a  and b: BTW your procedure doesn't use the parameter a)?

The sequence of binary digits you want is a binary de Bruijn sequence of order n.  Following
<http://en.wikipedia.org/wiki/De_Bruijn_sequence>, here's a construction using the GraphTheory package in the case n=5.

> with(GraphTheory);
    V:= [seq(seq(seq(seq(cat("V",a,b,c,d),a=0..1),b=0..1),c=0..1),d=0..1),e0,e1]:
# I'm using the two extra vertices e0 and e1 to take care of the fact that GraphTheory does not allow loops as edges
    E:= remove(t -> (t[1]=t[2]),{seq(seq(seq(seq(seq([cat("V",a,b,c,d),cat("V",b,c,d,e)],a=0..1),b=0..1),c=0..1),d=0..1),e=0..1)})
 union {["V0000",e0],[e0,"V0000"],["V1111",e1],[e1,"V1111"]}:
    G:= Graph(V,E):
    IsEulerian(G,'T');
    T:= subs(e0=NULL,e1=NULL,T);
    S:= map(parse,[seq(op(1,T)[i],i=2..5),seq(op(i,T)[5],i=2..nops(T))]);

S := [0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0]

Perhaps you mean something like ?examples,binarytree

> solve({cos(x)*sin(y)=0, y=2*sin(x)},AllSolutions);

,

That square-root term makes a big difference: in particular, it makes your equation nonlinear. 
Actually, you can get solutions if you express your pde in a way that eliminates the square root:

> pde:= (-alpha^2*f(x,y,z)-diff(diff(f(x,y,z),x),x)-diff(diff(f(x,y,z),y),y)-diff(diff(f(x,y,z),z),z))^2 = a^2*(diff(f(x,y,z),x)^2+diff(f(x,y,z),y)^2+diff(f(x,y,z),z)^2)

> pdsolve(pde);

I don't understand what you're trying to do.  You have one equation, with x the only variable in it (and in particular no w).  That equation has two positive real solutions, x=1 and x=2.056295775.  Once you have that other solution, why are you trying fsolve again (and without the "avoid" option) for each value of w, and why are you surprised that it gives you 1 each time?  If you want the other solution for x, just save that value in a variable and use it, e.g.

> X:= fsolve(equ1, x, avoid = {x=1});
   {seq({w = W, x = X}, W = wvals)};

Or did you have something else in mind?  Perhaps a different equation that has a w in it?

The color option, it seems, only affects the symbols plotted, not the lines joining them.  You could do this, however
(with de as in acer's answer):

> InitialValueProblem(de, t = 1, output = plot, numsteps = 10, method = euler, comparewith = [[rungekutta, rk4]]);

To show an antiderivative with +C, you have to add the +C yourself. But for the other problems, you can use definite integrals.  Alternatively, you could use dsolve.  For example:

> dsolve({diff(x(t), t) = v(t), diff(v(t), t) = cos(t), x(0) = 3, v(0) = 0});

M:= <<a,d>|<b,e>|<c,f>>;
Matrix(2,3,combinat[randperm](convert(M,list)));

As you say, this is an oscillatory integral, and it's not unexpected that this should cause problems for numerical integration techniques.  It's almost the same as the first question in the SIAM 100 Dollar - 100 Digits Challenge, see
http://en.wikipedia.org/wiki/Hundred-dollar,_Hundred-digit_Challenge_problems and http://www.math.ubc.ca/~israel/challenge/challenge.mws (but missing the factor 1/x, which makes it easier).