You operate with matrices whose elements are also matrices. When you multiply these objects, Maple does not understand that the individual elements of these objects should also be multiplied as matrices. Workaroud is that probably you should write this multiplication explicitly

See an example and a workaround:

A:=Matrix(2,{(1,1)=Matrix([[1,2],[3,0]]),(1,2)=Matrix([[-1,0],[1,3]]),(2,1)=Matrix([[1,2],[4,0]]),(2,2)=Matrix([[-1,0],[1,5]])});

B:=A^2;  # Maple does not know what to do next

B1:=Matrix(2, (i,j)->add(A[i,p].A[p,j],p=1..2));  # A workaround

restart;

Exponents1:= proc(p, x::name)

 local L, k, i, j;

k:=0:

for i in op(p) do

for j in op(i) do

if has(j,x) then k:=k+1; if nops([op(j)])=1 then L[k]:=1 else L[k]:=op(2,j) fi; fi;

od; od;    

convert(L, list);

end proc:

Your example:

f:=5*A^4*B^3*C^7+3*A^2*B^1*C^7+31*A^3*B^6*C^11;

Exponents1(f,A);       

Exponents1(f,B);

Exponents1(f,C);       

 Addition.

Unfortunately Maple  (Maple 2015.1) does not preserve the order of terms in the polynomial. If you want the order  of the exponents has been saved, you can specify the polynomial as a list of its terms, as follows

f:=[5*A^4*B^3*C^7, 3*A^2*B^1*C^7, 31*A^3*B^6*C^11];

Procedure code remains unchanged.

Another addition. 

I read the acer's solution and noticed that my solution is not working properly, if the polynomial consists of a single term. I did not bother to correct my solution as acer's solution easier and more understandable.

restart;

SplitScan1:= proc(f, L::list)

local M, m;

M:=[seq(`if`(f(L[i],L[i+1]), i, NULL), i=1..nops(L)-1)];  # List of positions of L in which it must be split

m:=nops(M);        

[L[..M[1]], seq(L[M[k]+1..M[k+1]], k=1..m-1), L[M[m]+1..]];

end proc:

The procedure is easy to modify if you want to find not the whole split list L, but some of its sub-lists, because the all information about the positions of split is in the pre-calculated list M. This can save a lot of time if the list  L  is long.

I fixed a few syntax errors in your code. 

ode := diff(y(x), x, x, x)+2*(diff(y(x), x, x))/x+2*(diff(y(x), x))/x^2-(25/8)*x^2*(7*x^10+52*x^5+16)*y(x)^7 = 0;

sol := dsolve({ode, y(0) = 1, (D(y))(0) = 0, ((D@@2)(y))(0) = 0}, numeric);

plots[odeplot](sol, [x, y(x)], x = 0 .. 1);

If you build a plot of a finite list of points in  style=point , then the points with the same values of the function can be identified by one color.

Simple example:

L := [seq([(1/6)*Pi*n, abs(sin((1/6)*Pi*n))], n = 0 .. 12)]:

M := ListTools[Categorize]((x, y) -> x[2] = y[2], L):  # Selection of points with the same ordinates

plot([M], style = point, symbol = solidcircle, color = [blue, green, yellow, red], scaling = constrained);

Another way - just use the formula for solution of  the equation  A1.A2=A3 . It does not require a call of  LinearAlgebra  package:

A2 = A1^(-1).A3;

The result is the same.

If you want to have a result without the binary variable  _B1  that simply substitute its values:

restart:

k:=combine((sin(x))^2);

sol:=solve(k=1/4, x, AllSolutions = true, explicit);

about(_B1);

sol1, sol2:=subs(_B1=0,%), subs(_B1=1,%);

 Addition: you can do the same with the other variable  _Z1 . Both replacements can be done simultaneously.

Done in Maple 2015.1

restart;

A := <a,b,c>:

a := 1: b := 2: c := 3:

B:=convert(A, list);

B;

                   B := [a, b, c]

                     [1, 2, 3]

The code was edited.

Use  Sum  instead of  sum  to prevent premature calculation.

restart;

Squaring:=proc(Expr)

local f, expr, E, n, r;

f:=op(0,Expr); expr:=op(1,Expr);

if f=Sum then

n:=op([2,1],Expr); r:=op([2,2],Expr);

return Sum(expr^2,n=r) + 2*Sum(Sum(subs(n=_m, expr)*expr, _m=op(1,r)..n-1), n=op(1,r)+1..op(2,r)) else Expr^2 fi;

end proc:

Example:

Squaring(Sum(sin(n*x), n=1..N));

value(%);

Another very simple example:

Squaring(Sum(a[k], k=1..5));

value(%);

Remark.  @Mac Dude  your term  sum(sum(E[m]*E[n],m=1..n-1),n=1..N)  is incorrect because should be  m<>n 

The code was edited.

restart;

seq([b,solve({x+b*y=0, b*x-y=10},{x,y})], b=10..20);

The use of strings can solve your problem. This code is easily rewritten as a procedure and used for similar tasks.

a:=2:  b:=4:  c:=4:  L:=["a","b","c"]:  L1:=map(convert, [a,b,c], string):

S:="a + b*c":  S1:=parse(S): T:=S:

for i to 3 do

T:=StringTools[SubstituteAll](T, L[i], L1[i]);

od:

convert(cat(S, " = ", T, " = ", S1), symbol);

A:=1.330169822*10^11/(6.552475529*10^9*a+7.554142204*10^9*b);

expand(``(numer(A)/10^9)/``(denom(A)/10^9));

For example:

Explore(eval(y*z+x, {x = x0, y = y0, z = z0}), parameters = [x0 = -10 .. 10, y0 = -10 .. 10, z0 = -10 .. 10]);

The list for f(x,y,z)  if  z=-10..10  (for example if x = 1, y = 2):

[seq(eval(y*z+x, {x = 1, y = 2}), z = -10 .. 10)];

@Rouben Rostamian  Your  z = 1/((1+x)*(1 + sqrt(x /a)*arccot(1/sqrt(a/x))))  does not coincide with the initial   z=1/((1+x)*(1+sqrt(x/a)*arccot(sqrt(a/x)))) . 

It seems Maple can not calculate this integral symbolically  even for concrete  a , but it is easily calculated numerically. This calculation is expedient to issue as a procedure with the parameter  a .

Integral:=a->evalf(Int(1/(1+x)/(1+sqrt(x/a)*arccot(sqrt(a/x))), x=0..infinity)):

We can use this procedure to evaluate the integral for the given numbers a and plot the integral vs a.

Examples:

Integral(0.1), Integral(1), Integral(5);

plot(Integral, 0..10, color=red, thickness=2);

If I understand your problem, here is the procedure to solve it. The last parameter of the procedure is N (by default N=100000) to prevent an infinite loop.

FirstReturn := proc(F::procedure, x::realcons, p::realcons, q::realcons, N::posint:=10^5)

local z, p0, q0, t, u, u0, z0;

z:=x; p0:=evalf(p); q0:=evalf(q);

for t from 1 to N do

u:=z; z:=F(z); u0:=evalf(u); z0:=evalf(z);

if (u0<p0 or u0>q0) and (z0>=p0 and z0<=q0) then return [z0, t]  fi;

od;

end proc:

Example of use. Let's try to solve the following problem. Lengthwise of the unit circle from the point (1,0) we shift by 1 (moving counterclockwise). Find the step in which we will be again in range 0 .. 0.0001 (from the point (1,0)).

FirstReturn(x->x+1-floor((x+1)/(2*Pi))*2*Pi, 0, 0, 0.0001);

                                    [0.0000602, 710]

It is known that if the shift (in this example it is 1) incommensurable with the length of the circle (in the example it is 2*Pi), the resulting set of points is dense on the circle.