There is no need to write the procedure choosing. Your example can be implemented in Maple as

L:=`$`~([w,b,g],4);
combinat:-choose(L, 4);

                L := [w, w, w, w, b, b, b, b, g, g, g, g]
 [[b, b, b, b], [b, b, b, g], [b, b, b, w], [b, b, g, g], [b, b, g, w], [b, b, w, w], [b, g, g, g], [b, g, g, w], [b, g, w, w], [b, w, w, w], [g, g, g, g], [g, g, g, w], [g, g, w, w], [g, w, w, w], [w, w, w, w]]

To solve the second problem, use the  ListTools:-Occurrences  command, for example:

e:=[b, b, b, w]: L0:=[w,b,g]:
[seq(s=ListTools:-Occurrences(s,e), s=L0)];

                             [w = 1, b = 3, g = 0]

Edit.

The procedure  Ellipse  does what you want:

restart;
Ellipse:=proc(a,b,phi)
local Vdot, dist, ell, F1, F2, M, MF2, MF1, varphi, c, m, t0, ELL, tp, po;
uses plots;
c := sqrt(a^2-b^2);
m:=2*arctan(c/b);
print(`maxϕ`=m);
if is(phi>m) then error "Should be phi<=`max&varphi;`" fi;
Vdot := proc(U, V) add(U[i]*V[i], i = 1 .. 2) end proc:
dist := proc(M, N) sqrt(Vdot(M-N, M-N)) end proc:
ell := x^2/a^2+y^2/b^2 = 1: 
c := sqrt(a^2-b^2):
F1 := [c, 0]: F2 := [-c, 0]:
M := [a*cos(t), b*sin(t)]:
MF2 := dist(M, F2);
MF1 := dist(M, F1);
varphi := arccos((MF1^2+MF2^2-4*c^2)/(2*MF1*MF2));
t0:=fsolve(varphi=phi, t=0..Pi/2);
print(parse("M")=eval(M,t=t0));
print(parse("MF2")=evalf(eval(MF2,t=t0)));
print(parse("MF1")=evalf(eval(MF1,t=t0)));
ELL := implicitplot(ell, x = -a-2 .. a+2, y = -2-b .. b+2, color = blue):
tp := textplot([[F1[], "F1"],[F2[],"F2"],[eval(M,t=t0)[],"M"]], 'align' = 'above', font=[Times,15]):
po := plot([F1,F2, eval(M,t=t0)], style = point, symbolsize = 12, symbol = solidcircle, color = red):
display([ELL, tp, po], scaling = constrained, size=[800,400]);
end proc:

Example of use:

Ellipse(7,5,Pi/4);

Edit. The procedure  Ellipse  text has been updated. I have added the calculation of the maximum angle  `maxϕ`  and an error message if  phi > `maxϕ`

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