If I understood your question correctly, then you want to learn how to build such graphs in Maple. Your graphs are usual broken lines. The simplest way to plot them in Maple is simply to specify a list of vertices of the broken line. Below, as an example, in Maple, the first graph is built:

restart;
L:=[[0.5,0.43],[0.525,0.32],[0.55,0.28],[0.575,0.52]]:
plot(L, thickness=2, size=[500,300], title="Wnf", labels=["b",``], font=[times,16], labelfont=[times,18]);

restart;
Digits:=20:
L:=[x^2/5-y^2/2=1,y=x^3-7/2*x^2+2*x]:
f:=solve(L[1],x);
g:=solve(L[2],x);
Sol:=evalf(RealDomain:-solve(L));
plots:-implicitplot(L, x=-2..4,y=-3..3, color=[red,blue], gridrefine=3);
Area=simplify(fnormal(int(g[1]-f[1], y=op([2,2],Sol[2])..op([2,2],Sol[1]))), zero);

                                Area =0 .75833670466253534952

restart;
for i from 1 to 2 do
for j from 1 to 2 do
T[cat(i,j)]:=2*mu*varepsilon[cat(i,j)]-add(2*mu/3*varepsilon[cat(r,r)]*delta[cat(i,j)],r=1..3);
print(S[cat(i,j)]=T[cat(i,j)]);
od: od:

Or (all indices in roman):

restart;
f:=(x,y)->parse(convert(cat(x,y),string)):
for i from 1 to 2 do
for j from 1 to 2 do
T[f(i,j)]:=2*mu*varepsilon[f(i,j)]-add(2*mu/3*varepsilon[f(r,r)]*delta[f(i,j)],r=1..3);
print(S[f(i,j)]=T[f(i,j)]);
od: od:

restart;
ode:=x = (diff(y(x),x)^2+1)^(1/2)*diff(y(x),x);
mysol:=y(x)=_C1+int(  sqrt(-2+2*sqrt(4*a^2+1))/2,a=0..x);
check:=odetest(mysol,ode);
simplify(combine(check)) assuming x>0;

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