When you rotate your 3D graph with the mouse, you can see the current values  of theta, phi and psi at the left of the context panel at any moment.
If you want to plot your 3D graph for specific values  theta, phi and psi, then use the  orientation  option (in degrees).

You are confusing simplifying an expression with solving an equation.
Here's a way "to simplify" your desired:

restart;
# Example 1
ode:=y(x)=diff(y(x),x)^3*y(x)^2+2*x*diff(y(x),x);
new_ode:=PDEtools:-dchange({y(x)=sqrt(u(x))},ode,{u});
new_ode1:=subs(diff(u(x), x)=freeze(diff(u(x), x)),new_ode);
u(x)=thaw(solve(new_ode1, u(x)));
``;
# Example 2
expr:=sqrt(u)=1/sqrt(u);
u=solve(expr,u);

I was able to open your file. Here is an example of plotting 8 graphs for 8  q-values in a for-loop:

restart; 
m := 0.1: c := 0.06667:  h := 0.1: 
eq1 := diff(u(t), t) = u(t)*((1-u(t))*(u(t)-m)-q*v(t)); 
eq2 := diff(v(t), t) = h*v(t)*(u(t)-v(t)+c)/(u(t)+c); 
dsys3 := {eq1, eq2, u(0) = 1/2*(1+m-q-sqrt((1+m-q)^2-4*(c*q+m))), v(0) = 1/2*(1+m-q-sqrt((1+m-q)^2-4*(c*q+m)))+c};

for  q from 0.3 by 0.01 to 0.37 do
dsol5[q] := dsolve(dsys3, numeric);
print(plots:-odeplot(dsol5[q], [t,u(t)-q], t=0 .. 10,  caption = typeset("q","=",q), size=[1000,400]));
od:

Obviously the right side of your equation  eq1  is  > 0  for any  r ,  i.e.  u(r)>0 . But then it is also obvious that the integral diverges, i.e.  the right side = oo  for any  r>0  and   any  u(r)>0 .

Just consider 2 cases:

Check:=simplify(odetest(sol,ode));
s:=csgn(ln(x+1)-2);
eval(Check,s=1); # not OK
eval(Check,s=-1); # OK
solve(s=-1) assuming x+1>0;

 
Addition. I noticed that adding the  assuming x+1>0  solves the problem directly:

restart;
ode:=(x+1)*diff(y(x),x)+y(x)^(1/2) = 0;
ic:=y(0) = 1;
sol:=dsolve([ode,ic],y(x));
Check:=odetest(sol,ode);
solve(Check) assuming x+1>0;

Edit.        

In fact, that is 5 and not 4.5. Look carefully. I changed a little tickmarks option:

restart;
R:=10:
r:=5:
plot3d( [ ( R+r(u,v)*cos(v))*sin(u),
          ( R+r(u,v)*cos(v))*cos(u),
            r*sin(v)
        ],
            u=0..2*Pi,
            v=0..2*Pi,style=patchnogrid,
            scaling=constrained,
          scaling=constrained,
coordinateview=[-15..15, -15..15,-5..5],
lightmodel=light3, viewpoint = circleleft, tickmarks=[default,default,10]);

You have not specified the version of your Maple. Perhaps in the latest versions everything works as you would like. In my Maple 2018, everything is the same as yours.
I don't see any other way how to code the desired legend using the  plot, plots:-textplot, plottools:-rectangle  commands. Of course, this is quite troublesome.
Look at an example below:

restart;
interface(version);
``;
p1:=plot(x^2,x=-2..2,color=red, style=pointline,symbol=asterisk,symbolsize=14,adaptive=false,numpoints=20):
p2:=plot([[1,-0.9],[1.4,-0.9]],color=red):
p3:=plots:-textplot([1.6,-0.8,x^2],font=[times,16]):
p4:=plot([[1.2,-0.9]],style=point,color=red,symbol=asterisk, symbolsize=14):
p5:=plottools:-rectangle([0.9,-0.6],[1.7,-1.1], style=line):
plots:-display(p1,p2,p3,p4,p5, size=[500,500]);

restart;
R:=Uo/sqrt(Go*Ho)=Fr;
subs(Uo=solve(R,Uo), Go*Ho/Uo^2);

The equality  Uo=solve(R,Uo)  is equivalent to the equality  R .

Here is a more direct and shorter solution. I used multiple assignment for all constants, removed indices  m  and  n  that were not used in any way, specified the matrix  S  and vector  V  of the system in a shorter way, and used the  LinearAlgebra:-LinearSolve  command. The result is returned as a vector  Sol , the entries of which we can give any names, for example  A  and  B :

restart;
h, C, rho, omega := 0.03, 1, 7800, 222:
S := <rho*h*omega^2+100,200; 280,rho*h*omega^2+444>;
V := -C*<555, 6665>;
Sol := LinearAlgebra:-LinearSolve(S, V);
A, B := seq(Sol);   

restart;
L:=proc(P,v)
if P=x then  x*y elif 
P=y then  1 elif
P=z then  z^2 elif
nops(P)=1 and not (P in v) then  0 elif
type(P,`+`) then  map(L,P,v) elif
type(P,`*`) or type(P,`^`) then if type(P,`^`) then  op(2,P)*op(1,P)^(op(2,P)-1)*L(op(1,P),v)  else  `+`(seq(`*`(({op(P)} minus {p})[])*`if`(p::numeric,0,`if`(p in v,L(p,v),op(2,p)*op(1,p)^(op(2,p)-1)*L(op(1,p),v))),p={op(P)})) fi; fi;
end proc:

# Examples
L(x*y + y*z,[x,y,z]);
L(x*y*z,[x,y,z]);
L(1/2*x*y^2+3*y^3*z, [x,y,z]);                               

More complex examples can be easily generated using the  randpoly  command, for example:

P:=randpoly([x,y,z], dense, degree=3);
L(P,[x,y,z]);

Edit.

See Help on  inttrans:-laplace  command and  dsolve,inttrans  help page. 

Your function should be defined as  piecewise  function:

f:=t->piecewise(t>=0 and t<5,exp(-2*t),t>=5,-t);
inttrans:-laplace(f(t), t, s);

restart;
# Define i1(t) and i2(t) as symbolic variable

# Given differential equation
ode1 := diff(i1(t),t) = 0.5*i1(t)-3*i2(t)+5*exp(-2*t);
ode2 := diff(i2(t),t) = 2*i1(t) - 6*i2(t);
odes := {ode1, ode2};

# Define initial conditions
cond1 := i1(0) = 1;
cond2 := i2(0) = -1;
conds := {cond1, cond2};

# Solution of system of differential equations
Sol := dsolve(odes union conds);

plot(eval([i1(t), i2(t)], Sol), t=0..5, color=[red,blue], legend=["i1(t)","i2(t)"]);

To plot one or more points, the simplest way will be to use the syntax  plot([A, B, ...], style=point) , where A, B, ...  are points given by lists of their coordinates. To specify the color, size and shape of symbols representing points, use the options  color =... ,  size = ...  and  symbol = ...  . To label points we use  plots:-textplot  command. Note that the function  y=x^2  is defined as a procedure  f:=x->x^2 . This is convenient when calculating the values of both the function itself and its derivative at specific points.

Below, as an example, we solve the problem:

Given the curve  y=x^2  and 2 points  A(-1,1) ,  B(2,4)  .
1. Check that these points lie on this curve.
2. Find the point on the curve  y=x^2 and the equation of the tangent line to this curve (at this point)  such that the tangent is parallel to the chord  AB .
3. Draw all the objects in one plot.

Solution:

restart;
f:=x->x^2:  A:=[-1,1]:  B:=[2,4]:

# Calculations
#
x1, x2, y1, y2 := A[1], B[1], A[2], B[2];
is(f(x1)=y1), is(f(x2)=y2); # Checking that points A and B lie on the curve y=x^2
Eq := D(f)(x0)=(f(x2)-f(x1))/(x2-x1);  # Mean value theorem
x0:=solve(Eq, x0);
y0:=f(x0);
k:=D(f)(x0); # The slope of the tangent
Tangent:=y=y0+k*(x-x0); # The equation of the tangent 
C:=[x0,y0];  # The point of tangency

# Plotting
#
Lines:=plot([f(x), [A,B], rhs(Tangent)], x=x1-1..x2+0.5, color=[blue,green,red]): # The plot of all the lines
Points:=plot([A,B,C], style=point, color=blue, symbol=solidcircle, symbolsize=12): # The plot of all the points
Labels:=plots:-textplot([[A[],"A"],[B[],"B"],[C[],"C"],[-1,3.5,y=f(x)],[-1,-1.5,Tangent]], font=[times,16], align={left,above}): # All the labels
plots:-display(Lines,Points,Labels, scaling=constrained, view=[x1-1..x2+0.5,-2..5], size=[500,500]);

When you apply thу  simplify  command to some equality, it is separately applied to the left and right sides, i.e.  simplify(A=B)  is equivalent to  simplify(A)=simplify(B)  (that's the design). 
To check some equality (an identity) use the  is  command instead of simplify:

restart;
simplify((a^2-b^2)/(a-b)=a+b);
is((a^2-b^2)/(a-b)=a+b);

If you really want to simplify your expression  f  to the expression  g , there is a workaround. We make simple substitutions, after which the numerator and denominator will be polynomials, then simplify. I wonder if the Mathematica can directly simplify  f  to  g  without such tricks.

restart;
f:=(9*(x^(-2/3*a))^2*exp(6/a*(x^(-2/3*a))^(1/2))^2*_C0^2-6*(x^(-2/3*a))^(3/2)*exp(
6/a*(x^(-2/3*a))^(1/2))^2*_C0^2*a+x^(-2/3*a)*exp(6/a*(x^(-2/3*a))^(1/2))^2*_C0^
2*a^2+18*(x^(-2/3*a))^2*exp(6/a*(x^(-2/3*a))^(1/2))*_C0-2*x^(-2/3*a)*exp(6/a*(x
^(-2/3*a))^(1/2))*_C0*a^2+6*(x^(-2/3*a))^(3/2)*a+x^(-2/3*a)*a^2+9*(x^(-2/3*a))^
2)/(3*_C0*exp(6/a*(x^(-2/3*a))^(1/2))*(x^(-2/3*a))^(1/2)-exp(6/a*(x^(-2/3*a))^(
1/2))*_C0*a+3*(x^(-2/3*a))^(1/2)+a)^2:
g:=x^(-2/3*a):
Subs:=[x^(-2*a*(1/3))=A^2,exp(6*sqrt(x^(-2*a*(1/3)))/a)=B];
f1:=subs(Subs, f);
f2:=simplify(f1) assuming A>0;
subs(A=sqrt(x^(-2*a*(1/3))), f2);
is(%=g);

restart;
A:=Matrix(3, [$1..9]);
L:=[indices(A, pairs)];
L1:=(rhs=lhs)~(L);
T:=table(L1);
T[7];