Another opinion is to use Statistics:-Tally: 

convert~(Statistics:-Tally(L), list);
 = 
       

If you just want future input to be in document mode by default (so existing cells remain unchanged), you can simply open the file in a text editor (like Notepad) and change the fifth line and the sixth line from 

 to 

Besides, the remaining options may be deleted, then Maple will use the default settings for this file. As an example, you can also change them to: 

By the way, to further convert existing content from worksheet mode into document mode, one has to manually edit some other code. For instance, the following text 

should be converted into 

Note that these are undocumented, so there is no guarantee that they will always work …. 
However, it is relatively easy to convert an existing document in document mode to worksheet mode in Maple: just save it in the Maple Classic Worksheet (.mws) format. 

Here is a somewhat inefficient approach: 

e1 := a/sqrt(tan(x+c__1)^2+1):
(_ -> simplify(ifelse(type(_, atomic), _, map(thisproc, _)), trig))(e1);
 = 
                     

But the advantage is that there is no need to understand the structure of the expression beforehand.

There exist at least two approaches, 

e1 := a/sqrt(tan(x + c__1)**2 + 1):
convert(e1, exp, 'simplifier' = rcurry(simplify, trig, exp)@normal, 'evaluate');
 = 
                     

simplify(simplify(convert(e1, 'sincos'), constant, sqrt), trig);
 = 
                     

even though neither of them looks elegant …. (By the way, I believe it should be a BUG that the output of simplify(e1, trig); remains unchanged.) 

If you set interface('typesetting' = "standard"):, the eqn will look somewhat shorter: . 

"A& B& \\textrm{The fundamental matrix has } &C ":
StringTools:-RegSubs("(.*)\\\\textrm\\{(.*)\\}(.*)" = "\\1\\begin{minipage}{\\linewidth}\\textrm{\\2}\\end{minipage}\\3", `%`);
 = 
  "A& B& \begin{minipage}{\linewidth}\textrm{The fundamental 

     matrix has }\end{minipage} &C "


StringTools:-RegSubs("(.*)\\\\textrm\\{(.*)\\}(.*)" = "\\1\\begin{minipage}{\\linewidth}\\textrm{\\2}\\end{minipage}\\3", "\\textrm{ WILDCARD\n* }");
 = 
        "\begin{minipage}{\linewidth}\textrm{ WILDCARD

          * }\end{minipage}"

 Unfortunately, Maple's regular expression engine does not offer certain features, and therefore there exists a problem: 

StringTools:-RegMatch("(.*)\\\\textrm\\{(.*)\\}(.*)", " A&B&{\\textrm{The fundamental matrix has {}}}&C ", '_0', '_1', '_2', '_3'):
StringTools:-IsBalanced~([_ || (0 .. 3)], "{", "}");
 = 
                   [true, false, false, true]

In my view, a better approach is using Python's `regex` library (in Maple) instead. 

Although inefficient (much slower than the direct nested loops), this is still a good idea: 

s := (n::posint) -> 
   add(add(1 /~ 
      mul~(ListTools:-FlattenOnce(
        combinat:-permute~(
         combstruct['allstructs']('Partition'(l), 
          'size' = 3))))), l = 3 .. n):
Warning, (in s) `l` is implicitly declared local

seq(s(n), n = 1 .. 10);
 = 
                    5  17  49  967  801  4523  84095
           0, 0, 1, -, --, --, ---, ---, ----, -----
                    2  4   8   120  80   378   6048 

time[real](s(100));
 = 
                             6.046

Why not just 

[LinearAlgebra:-Column](Matrix(3,4,'symbol'=m),[1..-1]);
         


whattype~(`%`);
[Vector[column], Vector[column], Vector[column], Vector[column]]

I'm not sure why, but there're at least 2 workarounds: 

RealDomain:-solve('identity'(eq,x),trial_solution_constants);

PDETools:-Solve(eq,trial_solution_constants,'independentof'=x);

I am unclear about why the latter command is so inefficient. (Actually, Mma carried out the same thing without any manual interference within 2 seconds!) 
Anyway, as regards your example, certain manual interventions appear necessary. 

restart;
B_poly := -(x + 4)*(x + 3)*(x + 2)*(x + 1)*(x - 1)*(x - 2)*(x - 3)*(x - 4):
g := B_poly/100000*(((B_poly - 669)/670)^136 - ((B_poly + 669)/670)^136):
f := -4347225/87808*x^8 - 17375/392*x^7 + 629491375/395136*x^6 + 266375/252*x^5 - 200677775/12544*x^4 - 3174625/504*x^3 + 11067842125/197568*x^2 - 53625/98*x - 126496075/4116:

Firstly, we can see that we only need test four (closed) intervals. 

rts := applyrule('RealRange(a::anything, b::anything) = a .. b', [RealDomain:-solve](B_poly >= 0, x));
 = 
              [-4 .. -3, -2 .. -1, 1 .. 2, 3 .. 4]

Next, we shall check whether g - f changes signs on these four intervals. (In theory, using sturm(sturmseq(g - f, x), x, a, b) should be enough; unfortunately, this built-in function seems to be less efficient as well (and so is fsolve(g - f, x, a .. b, 'maxsols' = 1)).) Here is a workaround. 

CodeTools:-Usage([seq](traperror(RootFinding:-RefineRoot(k, expand(g - f), x)), k = rts)); # treated as open intervals 
 = 
memory used=11.05MiB, alloc change=0 bytes, cpu time=515.00ms, real time=520.00ms, gc time=0ns

        ["given interval does not contain any roots", 

          "given interval does not contain any roots", 

          "given interval does not contain any roots", 

          "given interval does not contain any roots"]

Now we know g - f is either positive or negative in each of (-4, -3), (-2, -1), (1, 2), and (3, 4), so we just need check if  at the endpoints and some interior point of each interval. However, since the inequality is not strict, in accordance with the continuity, it is enough to consider one interior point of each interval. 

andseq(eval(g - f >= 0, x = stats['describe', 'mean'](convert(k, list))), k = rts); # verify at the midpoint 
 = 
                             false

Accordingly, the region {B_poly ≥ 0, g - f ≥ 0} is empty. 

These computations can be finished in one second on my low-end computer. 

First and foremost, “:=” should be used instead: 

cont_real := Re(continuous_solution):
cont_imag := Im(continuous_solution):
disc_real := Re(discrete_solution):
disc_imag := Im(discrete_solution):

To display a list of 2-D points, “plots:-pointplot” should be used: 

p1 := plots:-pointplot(x_values, cont_real, 'color' = "blue", 'legend' = "Continuous - Real"):
p2 := plots:-pointplot(x_values, disc_real, 'color' = "green", 'legend' = "Discrete - Real"):
p3 := plots:-pointplot(x_values, cont_imag, 'color' = "red", 'legend' = "Continuous - Imag"):
p4 := plots:-pointplot(x_values, disc_imag, 'color' = "orange", 'legend' = "Discrete - Imag"):

However, since some data in `disc_real` and `disc_imag` are too large, it is better to draw "log10~(disc_real)" and "log10~(disc_imag)" instead.

How about

is_symbol_inside_func_only := proc(expr::anything, f::name, y::symbol, ` $`)::truefalse; 
	type(applyrule('conditional'(f(_x::anything), _depends(_x, y)) = `tools/gensym`(f), expr), freeof(y)) 
end: 

expr:=3*ln(1+y)+ln(3*y)*y+ln(y)+cos(7*y):
is_symbol_inside_func_only(expr,ln,y); #should return false

expr:=3*ln(1+y)+ln(3*y):
is_symbol_inside_func_only(expr,ln,y); #should return true

expr:=ln(y)+ln(3*y)+cos(y):
is_symbol_inside_func_only(expr,ln,y); #should return false

expr:=3+cos(y):
is_symbol_inside_func_only(expr,cos,y); #should return true

expr:=y+ln(y):
is_symbol_inside_func_only(expr,ln,y); #should return false
 = 
                             false

                              true

                             false

                              true

                             false

Sometimes it would be better if there exists a special cbrt function. However, here you can use surd(x, 3) directly. 
Edit. Note that 

showstat(RealDomain:-`^`, 3);

RealDomain:-`^` := proc()
       ...
   3   clean(`assuming`([proc( s, n )
           if type(n,'fraction') then
               surd(s,denom(n))^numer(n);
           else
               s^n;
           end if;
       end proc(_passed)],['real']))
end proc

So it is somewhat unnecessary to with(RealDomain, `^`):.

Alternatively, 

r1 := -1 <= x and x <= 0:
r2 := 0 <= x and x <= 1:
convert(solve(r1 or r2, {x}), `and`);
 = 
                       -1 <= x and x <= 1

Some help pages provide a rough indication of basic methods and algorithms used for a particular purpose inside. Sometimes setting infolevel[…] to a positive integer (or just using showstat(…)) is also conducive to seeing such information, but it may still be difficult to determine which formula was used simply by reading the source code.
For instance, according to 

showstat((combinat::Partitions)::NumberOfPartitions, 13 .. 14):

(combinat::Partitions):-NumberOfPartitions := proc(n, {method::identical(recursive,dp,hrr,auto) := 'auto'})
local pn;
       ...
  13     userinfo(1,':-numbpart',"Using Hardy-Ramanujan-Rademacher series algorithm");
  14     pn := combinat:-Partitions:-NumberOfPartitions:-numbpart1_hrr(n)
       ...
end proc

, Maple by default make use of the so-called Hardy–Ramanujan–Rademacher method for larger n as well. However, the artificial benchmark - partition function claims that