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

Yes, there is a way: 
 

Download Don't_Forget_Integrals!.mw 
Compare:

I wonder why by default Maple is unable to calculate the second antiderivative mentioned in the link given by the OP (Note that the substitution x = 3 - 1/(t - 1/8) can reduce some degree.) in terms of elementary functions. According to Integration Methods - method=pseudoelliptic, now Maple is capable of integrating it heuristically, but here Maple simply fails to do so (unless using method=Trager (yet the result is still less nice)).  (Besides, Mathematica is also rumored to discover human-friendly solutions to the so-called pseudo-elliptic integrals: https://github.com/stblake/algebraic_integration.)

restart;
_seed := 1234:
Warning, the use of _seed is deprecated.  Please consider using one of the alternatives listed on the _seed help page.

M := LinearAlgebra:-RandomMatrix(4, 'generator' = 0 .. 1., 'datatype' = float[4]);
st := time():
linalg:-exponential(M, -I*t):
time() - st;
                             0.141

st := time():
:-LinearAlgebra:-MatrixExponential(M, -I*t):
time() - st;
                             1.828

M := LinearAlgebra:-RandomMatrix(8, 'generator' = 0 .. 1., 'datatype' = float[4]);
st := time():
linalg:-exponential(M, -I*t):
time() - st;
                             4.219

st := time():
:-LinearAlgebra:-MatrixExponential(M, -I*t):
time() - st;
                             24.594

So the LinearAlgebra package does not always seem to be more powerful and efficient in doing linear algebra calculations?!

In theory, this should be a workaround: 

restart;
use MmaTranslator:-Mma in MapAll(Chop, 1.378834798932344e-15*I*t) end;
# Nevertheless, it does not work as expected. 
Error, (in MmaTranslator:-Mma:-MapAll) too many levels of recursion

But strangely, this works well: 

use MmaTranslator:-Mma in map(Chop, 1.378834798932344e-15*I*t) end use;
 = 
                               0

Though @acer's approach works here, the MmaTranslator:-Mma subpackage does not seem to allow us to go beyond the basics.

You may have to use: 

restart;
(proc(` $`)
	local `&*`;
	redefine(`&*`, ['associative']); # Regretfully, “setattribute(`&*`, 'associative');” does not work here.
	use `*` = `&*` in
		local expr := -(r0+Delta_r)^2*(46*r0-41*Delta_r)*r0^5;
		subsop(1=a,2=b,3=c,4=d, expr)
	end;
	eval(`%`, `&*` = `*`)
end)() = 
                            a b c d

The key is 'associative' (instead of `&*`). 

It is the exponential integral.

Another way may be 

DETools:-DEplot(diff(y(x), x) = abs(1/sqrt(1 - x^2)), y(x), x = -.99 .. .99, y = -1.6 .. 1.6, [y(0) = 0], 'arrows' = 'fish'):

According to the help page, “If the third argument n is present then it specifies the "truncation order" of the series calculations”, so you may simply use: 

series(BesselK(4, x), x = 0, 0);
                       (-4)      (-2)    / 0\
                   48 x     - 4 x     + O\x /

To avoid explicit "*" or "·", I think that you can use 

subs("⋅" = " ", InertForm:-Display((* an inert-form expression *), 'inert' = false));

As for the second point, maybe you can use 

Sol2 := InertForm:-Display((`%/`@`?[]`@[primpart*numer, denom]*primpart)( Sol2 ), 'inert'=false):

How about

(convert~)~(convert~(convert(x*y*z*z*w*w*w, list, `+`), list, `*`), list, `^`);
 = 
               [[[x, 1], [y, 1], [z, 2], [w, 3]]]

?

According to the documentation of Cmd.exe, if the directory path, files, or any information you supply contains spaces, you must use double quotation marks around the text. So a potential workaround is: 

restart; kernelopts(opaquemodules = false):
Compiler:-Build:-buildVarDB["X86_64_WINDOWS"]["LDLIBS"] := StringTools:-RegSubs(StringTools:-Join([convert(curry(sprintf, "(%s%s.*%s)"), function, StringTools:-Escape~([kernelopts('bindir'), kernelopts('dirsep'), ".lib"], 'regexp')) $ 3], " ") = "\"1\" \"2\" \"3\"", Compiler:-Build:-buildVarDB["X86_64_WINDOWS"]["LDLIBS"]()):
# Thereafter `Compiler:-Compile` will work. 

I am not sure if this works on your computer, but at least this works for me.
Note that on Windows platforms, the default external C++ compiler distributed with Maple is Clang, so one can download the latest release of the LLVM-based MinGW and then replace the subfolder $MAPLE_ROOT/llvm manually. 
 

Download 236724-How-Can-I-Get-The-Compiler-To-Work.mw