The command print never produces programmatically accessible output; the only way to access it is copy and paste.

I can't see any reason why you'd want your output as a set, although it's possible that you have a good reason.

Here are two syntax alternatives for the creation of a column Vector. The first two lines require relatively recent versions of Maple:

f:= x-> local y; fsolve(x*y=1):
R1:= <(for x from 1 by 0.2 to 2 do f(x) od)>;
R2:= <seq(f(x), x= 1..2, 0.2)>;
 

You need to correct the syntax in ode2. After that your dsolve with D(u)(0)=0 works perfectly:

ode2:= -v + mu*diff(r*diff(u(r), r), r)/r = 0:
dsolve({ode2, u(R)=0, D(u)(0)=0});
                                 2      2
                              v r    v R 
                       u(r) = ---- - ----
                              4 mu   4 mu

My corrections are shown in green. And, just to improve readability, I removed three levels of unnecessary parentheses.

It must be noted that while D(u)(0) = 0 implies u(0) is finite, the converse is far from true. D(u)(0) = 0 is a quite strong condition. 

The syntax for a multiple integral requires square brackets for the 2nd argument (the variables and ranges), like either of these:

simplify(Int(f(x,y,z)*Norm(n), [u= 1..2, v= 1..4]));
int(f(x,y,z)*Norm(n), [u= 1..2, v= 1..4], numeric);

However, your surface is undefined along the line y=3/4*x, which unfortunately passes through the rectangle (x= 1..2) x (y= 1..4). By some rearrangements of the variables, you can persuade int to give an answer, but that answer is Float(undefined).
 ​

This seems to me to be a bug in solve. The command eliminate, already shown by acer, is an algebraic generalization of solve, although it has far fewer options.

eqs:= {x*y^2-x, x*sin(Pi*y)}:
map2(eliminate, eqs, {x,y}); #Solve for each variable independently
         {[{x = 0}, {}], [{y = -1}, {}], [{y = 1}, {}]}

#Any return from eliminate whose 2nd element is not the empty set
#is not a complete solution (in the sense of solve). The possible
#presence of these incomplete solutions is what I mean by "algebraic
#generalization". There are no such solutions in this case; I include
#the next line only for the sake of completeness.
(C,InC):= selectremove(s-> op(2,s)={}, %);
      C, InC:= {[{x = 0}, {}], [{y = -1}, {}], [{y = 1}, {}]}, {}

#Remove the superfluous empty sets:
op~(..-2, C);
                  {{x = 0}, {y = -1}, {y = 1}}

#If you want to include the free variable in each solution:
map(solve, %, {x,y});
       {{x = 0, y = y}, {x = x, y = -1}, {x = x, y = 1}}

We must assume that discont returns all of the discontinuities, although it's pretty easy to have one's faith in this shaken. IMO computer algebra is not possible (at this time at least) without tacitly assuming that functions are piecewise continuous. 

The "flavors" used in RandomTools:-Generate can be nested arbitrarily. So, using the Vector flavor,

V:= RandomTools:-Generate('Vector'(variable(length = 1), 4));

The mul command will work for explicitly specified values of n. But to express the product for arbitrary n, you need to use product. Note how I iterate through the odd denominators only by using 2*a+1. Unlike (recent versions of) mul, the product command will not accept a "step" such as 2 as a 3rd argument. 

h:= x-> sin(x)/x:
J:= Int(product(h(x/(2*a+1)), a= 0..n), x= 0..infinity) 
    assuming n::nonnegint; 
Maple won't evaluate that integral; or at least it wouldn't do it in an amount of time that I was willing to wait. But it will do it for specified nonnegative integer values of n:

value(eval(J, n= 7));

For every n that I tried, it returned Pi/2. I suspect that there is a fairly easy proof that it's always Pi/2.

This seems too obvious---you must have considered it yourself: Is there any reason to not do simply

max(abs~(F-G))

where F is the matrix for f and G is the matrix for g?

Aposthropes (ASCII code 39) are not special in strings. The only special ASCII characters are double quotes (ASCII 92) and backslash (ASCII 34). These can also be included in strings; they simply need to be escaped with backslash if they're hard-coded.

The error message that you show is typical when using sum, and you mentioned sum in your title. So which did you actually use? The correct command is add, not sum. Furthermore, the eval and indexing are not needed. You should be able to replace the entire command with

evalDG(add(L));

Divide the integration interval symbolically, like this:

(h,f):= (x-p, x^2+p*x):  (a,b):= (0,10):
sol:= solve(h = f, x);
add(int~([f,h,f], x=~ [a..sol[2], sol[2]..sol[1], sol[1]..b]));
simplify(%);
                              (3/2)              
              1 / 2          \               1000
              - \p  - 6 p + 1/      + 50 p + ----
              6                               3  

This solution is still correct even if the solution points sol are not in the interval 0..10.

Any rtable (i.e., Array, Matrix, or Vector) can be shuffled with randperm without any need to directly access or refer to its entries or any need to directly refer to its dimensions with this short procedure:

Shuffle:= (R::rtable)->
    R[(combinat:-randperm@`[]`@`$`)~([rtable_dims](R))[]]
:

It looks to me as if it did name the function A. And what is there to optimize in returning a 2x2 matrix of zeros? You can just ignore those warnings.

@adel-00 You wrote:

• u := proc (x, t) options operator, arrow; tanh(2*Pi-t-x) end proc;
  v := proc (x, t) options operator, arrow; tanh(2*Pi+t+x) end proc;
  A := 1/4*(-2*u(x, t)*(diff(v(x, t), x))+2*beta1*u(x, t)*(diff(u(x, t), t, t))+beta2*u(x, t)^2*(u(x, t)^2+v(x, t)^2)-2*alpha1*u(x, t)*(diff(v(x, t), t, t, t))-alpha2*u(x, t)*(u(x, t)^2+v(x, t)^2)*(diff(v(x, t), t))+2*v(x, t)*(diff(u(x, t), x))+2*beta1*v(x, t)*(diff(v(x, t), t, t))+beta2*v(x, t)^2*(u(x, t)^2+v(x, t)^2)+2*alpha1*v(x, t)*(diff(u(x, t), t, t, t))+alpha2*v(x, t)*(u(x, t)^2+v(x, t)^2)*(diff(u(x, t), t)));
  f := (x,t) -> A:
  r1:=int(int(A,x=-10..10),t=-3..3,numeric):

Maple is attempting a symbolic integration on the inner integral because only the outer integral specifies numeric. This is what takes a long time. To appreciate the complexity of your A consider that expand(A) has 3730 terms (although I'm not saying that that is how the symbolic integrator is looking at it).

To integrate both integrals numerically, you need to factor out the symbolic constants. This is fairly easy to do:

Coeffs:= proc(e::algebraic, V::{list,set}(name))
local T, C:= coeffs(e, V, T);
    table([T]=~[C])
end proc
:
V:= [alpha1, alpha2, beta1, beta2]:
C:= Coeffs(A, V):
r1:= add(v*int(C[v], [x= -10..10, t= -3..3], numeric), v= [indices](C, nolist));
        r1 := -11.78614736 - 7.805610537*alpha2 + 0.6218523429*alpha1 + 
              100.4082421*beta2 - 7.650147471*beta1

This result only takes about a second.

By the way, your command f:= (x,t)-> A is unnecessary and also incorrect. If it were needed, the correct command would be f:= unapply(A, [x,t]).

[Edit: Code updated to accomodate terms with mulitple symbolic coefficients. Also, my earlier code wasn't reporting the constant term.]

I would ignore that help page statement about algsubs being a generalization. It is a highly specialized command that should only be used when it's specifically needed (which is rarely). In general, use subs or eval and do not even think of using algsubs. Also, subs and eval are far more efficient.