That is a very good question.

Some thoughts:

For the case you shared, I think, odetest cannot simplify the residual for some reason. I tried a few commands on the residual and then did something manual. The manipulation below involved a squaring step, which "removes" a negative sign

A+B=0 ->A=-B -> A^2=B^2 -> which is true/equal.

I interprete the loss of information by the squaring step in a way that the residual can be zero under certain conditions but not in general (depending on assumptions on a). Maybe that is the reason why odetest does not return zero (or there are limitations in algorithmic simplification).

I do not see such a loss of information in the derivation of the "root ode"s. Therefore the solutions of the "root ode"s should be valid for the original ode, because:
It should be possible to work backward from the "root ode's" to the original ode. Each algebraic step generates a new ode. At one point a squaring operation is required where the two branches of the "root ode"s collapse to a single ode (here the information about of the two branches is "lost", i.e. which solution belongs to which branch). If the  solutions of the "root ode"s satisfy the ode's before squaring, the ode after squaring should also be satified by the solutions. (I don't see a reason why squaring an equation compromises equaltiy.) Therefore the solutions of the "root ode"s should be valid for the original ode.

Download residual.mw

dsolve does not necessarily use in the first place what odeadvisor suggest.

dsolve goes through a lists of methods and stops as soon as a solution is found. In this case dsolve starts with

dsolve(ode,`dsolve/methods`[1, high_degree])

and stops when the method d'Alembert was successfull. The list containing the method "separable" was simply not executed. Maybe it was next in line.

dsolve(ode,`dsolve/methods`[1])

The returned result of "separable" seems to be correct but can be verified only by hand

odetest(sol_2,ode)=0;
map(exp,%);
simplify(%)

Why odeavisor returns only one classification is unclear

DEtools:-odeadvisor(ode,[dAlembert])
                          [_dAlembert]

Maybe its because odeadvisor suggest in this case only the method that Student:-ODEs:-ODESteps uses.

I remember a case you mentionned where two methods where suggested. On this case ODESteps failed.

Because the last two plots are animations. You have to run the animations

The document is not editable

To edit check the box.

I cannot give an answer why the global Typesettings are not used inside plots (acer probably can). I assume that the label option simply ignores this.

Here are two manual ways to get the y axis labeled as you want.
 

DEtools:-DEplot([eq1,eq2],[X1(x),X2(x)],x=0..100, 
            X2=-4..4,X1=-4..4, 
            axes=boxed,
            linecolor = red,        
            labels=[y(x),`#mrow(mi("y"),mo("'"),mfenced(mi("x")))`]);

DEtools:-DEplot([eq1,eq2],[X1(x),X2(x)],x=0..100, 
            X2=-4..4,X1=-4..4, 
            axes=boxed,
            linecolor = red,        
            labels=[y(x),typeset(y,"'(",x,")")]);

You get an indication in this (very interesting case) by giving solve more freedom (by not indicating variables) to search for other solutions or you have to use the "right" variables straight away.

Download solve-nonzero-solutions_with_solve.mw

Try to put the solve command into brakets. This way simplify gets all 3 solutions. 

(At the moment, I do not have Maple runnig to test this)

complexplot3d as you have used it cannot plot two surfaces (two complex expressions) descirbed with two parameters. You might wanted to do this

plots:-complexplot3d([Re(solnum[1]),Im(solnum[2])], x = -50.. 50, t = -50..50);

where the first expression represents the real part and the second the imaginary part.

Otherwise you can use one of

plot3d(map(Re, solnum), x = -50 .. 50, t = -50 .. 50);
plot3d(map(Im, solnum), x = -50 .. 50, t = -50 .. 50);
plot3d(map(abs, solnum), x = -50 .. 50, t = -50 .. 50);
plot3d(map(argument, solnum), x = -50 .. 50, t = -50 .. 50);

or you replace the parameters x and t by a complex parameter. For example

subs(x = Re(z), t = Im(z), solnum);
plots:-complexplot3d(%[1], z = -50 - 50*I .. 50 + 50*I);

Because it's not equal. Only if you square it.
 

is(sqrt(-2) = sqrt(2)*I);
is(sqrt(-2) = -I*sqrt(2));
                              true

                             false

One keyelement to understand the behaviour of Maple are double brackets ()(). This is very important but sparsely documented. You will probably not find it in the help system. You need ()() for function composition. See ?operators,functional and ?examples,functionaloperators.

The later explains differences between 1D and 2D.

(1.5) above not really explains why in your expression (4) the independed varaible x is dropped. I assume that sin^-1 is not interpreted as an operator/function and therefore the argument is dropped.

For more on ()() see acers answers here.

Concerning @@-1:

It makes Maples output for negative and postive exponents consistend.

I do not see the equaltiy you are claiming.
Assuming b real leads to roots for x that depend on beta

Only x=Pi/2 establishes equality, which makes B__0=-beta^5. 

The first call to TestDimension is ambigeous and allows many valid combination of dimensions. The second input does not (read the output backward) because x was replaced by 2.

From the help page:

In some situations, there is some ambiguity in how units are assigned, for example in the expression
                      "5*Unit('kg') + x*y"

if no extra information is given. In such a case, Maple will give the fully general value if you ask for
                     "output = dimensions"

This output is sometimes hard to understand, but it is explained in the text above.

I addition to acer answer:
In Units:-Standard Maple assumes that x is unitfree and therefore returns in the error message `1` for the unit of x.

In Unit:-Simple Maple does not assume anything (edit:) for the whole session for 5 * Unit(m) +x. The unit of x is left open and no dimensional checks are performed (edit: between inputs). unless x is assigned to a value.

Edit: Assigning x to a value with unit, does not allways trigger consistency checks
 

restart;with(Units:-Simple):
6*Unit('s') + x;
x:=2*Unit(m);
%%;
eval(%);
simplify(%);
                         6 Unit(s) + x

                         x := 2 Unit(m)

                     6 Unit(s) + 2 Unit(m)

                     6 Unit(s) + 2 Unit(m)

Error, (in Units:-Simple:-+) the following expressions imply incompatible dimensions: 6*Units:-Unit(s)+2*Units:-Unit(m)

What is currently not possible in Maple is explictily assuming that x is of a certain dimension. With that assumption Maple could check if operands of a sum are of the same dimension or collect a common factor from sums containing units.

To sure if this answers fully your question but if you in integrate the lefthand side of expression (2) over r and add an integration constant you will get an equivalent result as shown in DESol.

Download DESol_Question_reply.mw

The problem with the integral is that ln(V), where V is a physical quantity,  does not evaluate to numerical values.
Here is why Maple can only provide an answer if n and m are and the integration variable are assumed dimensionless:
 

Updated worksheet: Download Units_Int_reply_2.mw

Former worksheet: Download Units_Int_reply.mw

Assumptions on names having a unit such as

       assuming V::[m^3]

would be of great value for mixed expression containing names/symbols and units since Maple could detect conflicts and errors