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];

Eq:=a^(5/2)*s^2*(D(f))(eta)^2*sqrt(nu)/(R*sqrt(a)*sqrt(nu)+eta*nu) = a^(5/2)*s^2*(D(P(s, eta*sqrt(nu)/sqrt(a))))(eta)/sqrt(nu);
Term:=select(has,rhs(Eq),P);
solve(Eq, Term);

It's better to implement  this as a procedure so that  Digits  is set automatically:

restart;
Evalf:=proc(x::realcons,n::posint)
local m;
m:=ilog10(evalf(x));
Digits:=m+n+3;
printf(cat("%.",n,f),x);
end proc:

Examples of use:

Evalf(123456789/123,6);
``;
Evalf(123456789/123,15);
``;
Evalf(Pi,3);