@torabi I do not understand what you mean by "other option".

Maybe you mean a numerical solution? In this case, you also get a zero solution:

restart;
sol := dsolve(eval({diff(r*(diff(u(r), r)), r)/r = 0, u(R0) = U2, u(h) = U1}, [h = 1, R0 = 5, U1 = 0, U2 = 0]),  numeric);

plots:-odeplot(sol, [r,u(r)], r = 1 .. 5, color=red, thickness=3);
 

@brian bovril  It works for me in Maple 2015.2, Maple 2016.2 and Maple 2017.3 properly. Try to do restart and run again.

PS. I uploaded your worksheet. Everything works as well.

@gmzsvsclk 

restart:
P:=product((x-j),j=0..11);
Q:=a*x^6+b*x^5+c*x^4+d*x^3+e*x^2+f*x+g:
R:=j*x^4+k*x^3+l*x^2+m*x+n:
T:=Q^2-R:
[seq(coeff(T,x^p)=coeff(P,x^p), p=1..12), coeff(T,x,0)=coeff(P,x,0)
];
solve(%);

Verification:

assign(%[1]):
'Q'=Q;
'R'=R;
solve(Q^2-R);

Edit.

@asa12  In Maple 12 replace the line  
S:=[seq(subs(var=~p,obj2), p=P)];
by the line
S:=[seq(subs(zip(`=`,var,p),obj2), p=P)];

As for other matrix transformations (in particular those you are writing about), other procedures or an extension of  IsEquivalent  procedure are needed.

@asa12  I did not understand the meaning of your question.

Similar problems arise only in the classic interface. In standard one polygon's filling works correctly for any simply-connected polygons (convex or concave). I added only one line to your code (in Maple 2017.3):

plots[display](<l | lM | H>,scaling=constrained, color=gray);

@das1404  The file is so big because of the animation. Try not to open this file, but just copy the following commands as text and run them in your Maple. I think this should work in Maple 7 (there are no references to the font here). Run not all the commands at once, but one at a time.

restart;
Picture := proc (L, C, N::posint := 100, Boundary::list := [linestyle = 1]) 
local i, var, var1, var2, e, e1, e2, P, Q, h; 
global Border; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then P[i] := op(L[i]) else 
var := lhs(L[i, 2]); var1 := lhs(rhs(L[i, 2])); var2 := rhs(rhs(L[i, 2])); 
h := (var2-var1)/N; 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := seq(subs(var = var1+h*i, [e*cos(var), e*sin(var)]), i = 0 .. N) else P[i] := seq([var1+h*i, subs(var = var1+h*i, e)], i = 0 .. N) end if else 
e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; P[i] := seq(subs(var = var1+h*i, [e1, e2]), i = 0 .. N) end if end if end do; 
Q := [seq(P[i], i = 1 .. nops(L))]; 
Border := plottools[curve]([op(Q), Q[1]], op(Boundary)); 
[plottools[polygon](Q, C), Border]; 
end proc:

L1:=[[[0.79*cos(t),-0.13+0.8*sin(t)], t=-Pi/2..Pi/2], [[0,0.67],[-0.83,0.67],[-0.83,0.63],[-0.75,0.63]], [[-0.75+0.14*cos(t),0.49+0.14*sin(t)], t=Pi/2..0], [[-0.61,0.49],[-0.61,-0.76]], [[-0.75+0.14*cos(t),-0.75+0.14*sin(t)], t=0..-Pi/2], [[-0.75,-0.89],[-0.83,-0.89],[-0.83,-0.93],[0,-0.93]]]:
L2:=[[[-0.14+0.66*cos(t),-0.13+0.717*sin(t)], t=-1.94..1.94], [[-0.38,0.54],[-0.38,-0.8]]]:

LetterD:=plots[display](plot([0.62,[-0.75+0.14*cos(t),0.48+0.14*sin(t),t=0..Pi/2]], x=-0.83..-0.75, color=white,filled), Picture(L2, color=white, [color=blue, thickness=5]), Picture(L1, color=yellow, [color=blue, thickness=5]), scaling = constrained):
LetterD;

N:=60: # The number of the frames
R:=seq(plottools[rotate](LetterD, 2*Pi*k/N), k=0..N):
plots[display](R, insequence=true);

@das1404  You can see here an example of creating and rotating in 2D the letter D. It seems you did not see this my answer. I'm wondering whether this file   LetterD1.mws  will work in Maple 7?

@waseem   Change the last line to

plots:-display(Vector([seq(plot([seq(eval(lambda(F), Nb=j), j=[0.1,0.2,0.3])], delta2=0.02..0.1, legend=[Nb=0.1,Nb=0.2,Nb=0.3], title=('F'=F)), F=[10,20,30])]), size=[900,500]);

@ganelon  From this specification nothing will change: the solution is unique for any  y>0

@acer  How did you know about [Gr,Br] ?

@Carl Love  You're right.

@das1404   If I need to draw several graphs on one plot, then I specify the common options (axes, scaling, view) only at the end of  plots[display]  command. In this case there will be no contradiction between the options.

@Preben Alsholm 

We can simplify a little more by

EXPR2:=combine(EXPR1, ln) assuming h>0, R0>0;
length(EXPR2);

@alfarunner  There are no problems here. The coefficients before rho hat and phi hat are the coordinates of the vector v(t) in this basis (here phi=t):

restart;
f:=unapply(<r(t)*cos(t),r(t)*sin(t)>, t):
v:=unapply(diff~(f(t),t), t):
V:=<cos(t),sin(t); -sin(t),cos(t)>.v(t);
simplify(V);  # The coordinates in the new basis (rho hat and phi hat)

I have never worked in the Physics package.